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ABSTRACT 

Scott, Robert Kenmer. Ph.D., Purdue University, Decenibe 
1972. Thermo chemical Nonequilibrium in Atomic Hydrogen at 
Elevated Temperatures. Major Professor: Frank P. Incropera. 

A numerical study of the nonequilibrium flow of atomic 
hydrogen in a cascade arc was performed to obtain insight 
into the physics of the hydrogen cascade arc. A rigorous 
mathematical model of the flow problem was formulated, incor- 
porating the important nonequilibrium transport phenomena and 
atomic processes which occur in atomic hydrogen. Realistic 
boundary conditions, including consideration of the wall 
electrostatic sheath phenomenon, were included in the model. 
The governing equations of the asymptotic region of the 
cascade arc were obtained by writing conservation of mass 
and energy equations for the electron subgas, an energy 
conservation equation for heavy particles and an equation of 
state. Finite-difference operators for variable grid spacing 
were applied to the governing equations and the resulting 
system of strongly coupled, "stiff" equations were solved 
numerically by the Newton-Raphson method. 

Parametric solutions were obtained with arc current, 
pressure and tube radius as independent variables. It was 
observed that nonequipartition of kinetic energy between 
electrons and heavy particles was a relatively insignificant 
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effect. However, although chemical equilibrium existed near 
the centerline of the arc, significant departure from this 
condition occurred in the wall region. The radial location 
marking the onset of chemical nonequilibrium shifted toward 
the centerline as the current and pressure decreased and as 
tube radius increased. 

Heat transfer in the hydrogen cascade arc is mainly due 
to radiation and the heavy particle conduction mechanism. 

Heat transfer by electron conduction and diffusion and the 
diffusion of ionization energy was found to contribute less 
than 2% to the total for all parametric conditions. Radia- 
tion becomes a dominant mode of heat transfer in the high 
pressure hydrogen arc, with radiative transfer accounting for 
over 50% of the heat loss for a pressure of 10 atmospheres. 

Speculations for the molecular hydrogen arc indicate 
that, if the molecular species were considered, higher elec- 
tron number density and temperature and lower heavy particle 
temperature would be observed in the wall region. However, 
it is expected that bulk parameters, such as electric field 
and the total and radiative wall heat flux, would not be 
substantially different for molecular hydrogen. 
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CHAPTER 1 INTRODUCTION 
1.1 The Cascade Arc 

In the past, much interest has been shown by researchers 
in the phenomena associated with high temperature gases. 

This interest centered on as trophysical and geophysical 
problems, but more recently attention has been focused on 
man-made devices which are capable of heating gases to high 
temperatures. One such device, commonly referred to as the 
cascade arc (a gas flowing through an electrical discharge 
confined by the cooled wall of a tube) , has attracted the 
attention of a number of theoreticians and experimentalists 
(see Figure 1-1). This theoretical study is concerned with 
the flow of hydrogen in the cascade arc. 

For purposes of discussion, the flow field in the arc 
constrictor is divided into three general regions (Figure 1-1) : 
(1) the entrance region, (2) the asymptotic region and (3) the 
field free region. The entrance region is characterized by 
significant variation of flow properties in the axial direc- 
tion, the asymptotic (fully developed) region by its lack of 
such property variation and the field free region by the 
absence of an external electric field. 

Due to the aforementioned characteristics of the 
asymptotic region, the flow properties determined at one 
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z 




Figure 1-1. Schematic Diagram of Cascade Arc. 
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axial location for a given set of operating conditions are 
applicable to all other axial locations (of the same region) 
and thus possess a ”semi-universal'' nature. This implies 
that the flow properties of the asymptotic region are inde- 
pendent of the cathode and inlet geometry as well as axial 
location, a fact which makes meaningful comparisons between 
theory and experiment possible. Also, due to their simpli- 
fied form, the asymptotic region governing equations can be 
used in conjunction with experimental measurements of cascade 
arc flows to experimentally determine high temperature gas 
transport properties. Therefore, due to the practical 
importance of this region, the present study of hydrogen 
flow through the cascade arc has emphasized the determination 
of conditions in the asymptotic region. 

Many previous studies of the cascade arc have incorpor- 
ated the assumption that the gas flowing through the arc is 
in local thermodynamic equilibrium (LTE) . (The ramifications 
of this assumption are discussed in the next section.) More 
recent studies [1, 2] of high temperature gas flows indicate 
that the LTE assumption is not valid in regions of high 
temperature and concentration gradients such as might exist 
near the wall of a cylindrical cascade arc. Therefore, for 
this particular study, the LTE assumption is not utilized, 
thus permitting a critical comparison of the equilibrium and 
nonequilibrium theories of the cascade arc for hydrogen. 
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1 . 2 Equilibrium and Nonequilibrium in the Cascade Arc 

In order to give the reader a proper perspective, the 
concept of equilibrium and nonequilibrium in cascade arc 
flow is discussed in this section. The term ’’complete 
thermodynamic equilibrium" (CTE) is used to describe the 
state of a gaseous system which has been isolated from its 
environment for a sufficiently long period of time. This 
period of time must be long enough (theoretically infinite) 
to insure that the gas is in thermal, mechanical and chemical 
equilibrium, thus implying the absence of temperature, pres- 
sure and concentration gradients, respectively. Therefore, 
any system in complete thermodynamic equilibrium cannot 
possibly transfer heat, momentum or mass across the bound- 
aries of that system. Since virtually all gaseous systems 
of interest do not satisfy the GTE requirement, it would 
seem impossible to use the concepts of classical thermo- 
dynamics (such as pressure, temperature, internal energy, 
etc.) to solve practical flow problems. For this reason it 
has become necessary to develop the concept of local thermo- 
dynamic equilibrium (LTE) to establish a more reasonable set 
of criteria for applicability of thermodynamic concepts to 
flow systems. Therefore, if a system is in LTE, thermo- 
dynamic concepts can be used to mathematically describe it. 

Before discussing the definition of LTE for this prob- 
lem it would be beneficial to review some of the microscopic 
conditions which would exist in a high temperature gas 
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(T “ 10,000®K) in CTE. For the purpose of concreteness, it 
is assumed that the gas is composed primarily of electrons, 
atoms, molecules and singly-charged atomic and molecular 
ions. With the aid of statistical mechanics [3], the follow- 
ing facts can be demonstrated; 


1. The particle velocity distribution function of each 

chemical species is Maxwellian, and the species kinetic 
temperatures are given in terms of the root mean square 


y kT. = y C.2 
2 1 2 1 


speeds of these distributions 
all species kinetic temperatures are equal. 


In CTE, 


2. The populations of the bound electronic states of each 
neutral and ionic species are given by the Boltzmann 
distribution function evaluated at a common excitation 
temperature, T . 

3. The electron number density contribution of each 
atomic and molecular species (due to ionization) is 
given by the Saha equation for that particular species 
ionization reaction. The total electron number density 
is given by the sum of these individual contributions 
and the Saha equations are evaluated at a common 
ionization temperature, Tj. 

4. The neutral atom number density is given in terms of 
the neutral molecule number density by the law of mass 
action applied to the dissociation reaction and evalu- 
ated at some dissociation temperature, T^. 
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5. The radiation intensity within the gaseous system 
satisfies the Planck function at the radiation tempera- 
ture, Tj^. 

6. In CTE, there is a unique temperature, T, which 
describes all of the aforementioned phenomena. 

Therefore, in CTE, T = T^ = T^^ 

throughout the entire system. 

For the purposes of this problem, LTE is said to exist 
when the first four conditions plus the condition 
T = T^ = Tg^ “ satisfied at each point in the 

flow field. No requirements are placed on the nature of 
the radiation field. Therefore, in LTE, the thermodynamic 
temperature (and hence all temperature dependent properties) 
is allowed to vary throughout the flow field, thus permitting 
the molecular phenomena of heat, mass and momentum transfer 
to occur. 

The conditions necessary for LTE to exist also suggest 
several possible modes of nonequilibrium for the gaseous 
flow system. Taken individually, these nonequilibrium con- 
ditions are: 

1. Non-Maxwellian velocity distribution of one or more 
of the different species. 

2. Inequality of two or more of the species kinetic 
temperatures (nonequipartition of kinetic energy) . 
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3. Excitation nonequilibrium, or nonequilibrium popula- 
tion of atomic and molecular excited states, and 

4. Chemical nonequilibrium, or nonequilibrium number 
densities of the different chemical species (electrons, 
atoms, molecules, atomic and molecular ions). 

Processes which tend to promote or suppress the above 
nonequilibrium conditions are discussed below in their 
respective order. 

1. Strong electric fields may cause a "drift" motion of 
electrons with respect to the mass-averaged gas velocity. 
If this drift velocity is of the same order of magnitude 
as the mean electron thermal velocity, then the electron 
velocity distribution will be non-Maxwellian. Also, 
ionization and three-body recombination reactions 
selectively absorb and release energy in the high energy 
wings of the electron velocity distribution function, 
thus tending to promote a non-Maxwellian electron 
velocity distribution. This nonequilibrium tendency 

is retarded by elastic collisions among the particles 
of each chemical species. 

2. The electric field, by selectively imparting kinetic 
energy to the electrons, is partially responsible for 
nonequipartition between electrons and "heavy" particles. 
In addition, diffusion of the relatively mobile electrons 
away from the hot central core of the arc contributes to 
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nonequipartition in the outer periphery near the wall. 
Finally, the ionization (three-body recombination) and 
collisional excitation (or deexcitation) processes 
tend to absorb (or release) electron energy in those 
regions where they are dominant. In contrast, the 
primary process through which equipartition o£ energy 
is restored is the occurrence of elastic collisions 
between the electrons and the heavy particles. 

3. For an optically thin gas, the emission of radiant 
energy causes a depopulation of the excited states, 
since the photoexcitation restoration mechanism is not 
present. The collisional excitation process (when 
dominant) will serve to restore excitation equilibrium 
to the heavy particles in the gas. 

4. The existence of large concentration and electron 
kinetic temperature gradients in the presence of 
finite ionization, dissociation and recombination 
rates is primarily responsible for chemical nonequilib- 
rium. For example, diffusion of electrons from the 
hot core of the arc toward the cool wall will elevate 
the electron number density in the wall region above 
the equilibrium value for the local electron tempera- 
ture due to finite ionization and three -body recombina- 
tion rates. Therefore, restoration of chemical equi- 
librium can only be realized in the absence of large 
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concentration and species kinetic temperature gradi- 
ents or in the presence of very large reaction rates. 

Earlier in this section it was stated that, if a 
gaseous system is in LTE, classical thermodynamic concepts 
can be used in its mathematical description. However, 
since the LTE assumption is not utilized in this study, 
the use of any thermodynamic concepts needs to be justified 
for this particular nonequilibrium situation. 

As mentioned in Section 1.4, nonequipartition and 
chemical nonequilibrium are treated in this study. However, 
it is assumed that Maxwellian velocity distributions do 
exist for the individual species present. The important 
consequence of this assumption is that species kinetic 
temperatures can then be uniquely defined. Therefore, 
thermodynamic properties which are a function of the species 
kinetic temperatures (such as enthalpy and internal energy) 
are useful concepts, even in this nonequilibrium situation, 
and the governing equations can be derived by the same 
methods that are used to derive the governing equations for 
LTE flow. 


1 . 3 Previous Studies of Cascade Arc Flow 
It is the purpose of this section to familiarize the 
reader with the type of theoretical and experimental tech- 
niques used in the past to obtain the flow characteristics 
of cascade arcs. This is not intended to be even a cursory 
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review of the cascade arc literature, but simply an overview 
of some of the more popular methods of investigating flow 
phenomena in such devices. 

One of the earlier treatments consists of writing an 
energy balance for the gas in the asymptotic region in terms 
of an electric field (Ohmic) heating term, a Fourier conduc- 
tion term and a radiation heat loss term. The basic assump- 
tion in the resulting energy equation (Elenbaas-Heller equa- 
tion) is that the gas is in local thermodynamic equilibrium. 
Analytical solutions to this equation have involved its 
simplification through use of a heat flux potential and/or 
approximated coefficients. Numerical solutions have also 
been obtained utilizing the equation in its basic form. 

Stine and Watson [4] proposed a simplified method of 
analytically treating flow in the entrance region of a 
cascade arc. Assumptions used in their theory include LTE , 
constant mass flux throughout the constrictor, no radiation 
heat loss and gas properties which are linearly dependent 
upon enthalpy. Due to the limiting assumptions of the Stine- 
Watson model, it is capable of yielding only rough estimates 
of flow properties in the entrance region of the cascade arc. 

Although there are a variety of approximate methods 
used to analytically obtain the flow characteristics of 
cascade arcs , it is impossible to accurately treat the 
problem by non-numerical methods. For this reason, Bower 
[5] elected to derive a rigorous model of arc flow in a tube 
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and to solve it numerically with a minimum o£ simplifying 
assumptions. He used an implicit finite-difference, march- 
ing scheme to obtain the LIE flow characteristics for argon. 
Bower’s LIE program has been modified for hydrogen flow by 
Greene [6] and will be used to compare results obtained in 
this study. 

As mentioned in the first section, Incropera and Viegas 
[1] investigated the existence and nature of nonequilibrium 
in an arc by means of a time-scale study and demonstrated the 
probable existence of thermal and chemical nonequilibrium in 
argon. They then proposed a non-LTE fluid flow model, 
similar to those which had been used by other researchers 
to compute high temperature, non-LTE gas flows in various 
devices and geometries. For the sake of background informa- 
tion, some of these non-LTE studies are discussed. 

Okuno and Park [7] investigated nonequilibrium, stagna- 
tion point flow of nitrogen over a hemispherical body. After 
transforming variables, they obtained a system of coupled, 
nonlinear ordinary differential equations which were not 
solvable by the usual finite-difference techniques due to 
strong coupling of the equations. They, therefore, utilized 
a "shooting” method in which certain boundary conditions 
were estimated, a marching solution was then effected from 
that boundary, and the calculated results were compared with 
conditions imposed at the other boundary. The method of 
solution then consisted of a "conversational" iteration with 
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the computer until convergence upon the correct set of 
boundary conditions was obtained. 

Kruger [8] utilized a flow model similar to that pro- 
posed by Incropera and Viegas in order to establish the 
accuracy of spectroscopic data obtained from argon and helium 
confined arcs . Rather than solve the governing equations 
directly, Kruger integrated them to obtain explicit expres- 
sions for electron number density, electron and heavy parti- 
cle temperatures , and the heat flux potential in terms of 
the measured values of electron temperature and number dens- 
ity. He found that measured and calculated values of elec- 
tron temperature and density were in good agreement, heat 
flux potentials were in fair agreement, and that thermal 
and chemical nonequilibrium do exist in argon and helium 
arcs. His conclusion is that ”a two- temperature ambipolar 
diffusion model yields satisfactory interpretation of the 
nonequilibrium behavior of confined arcs in atmospheric 
pressure argon and helium.” 

The most rigorous theoretical prediction of cascade 
arc flow known to the author was completed by Clark [9] in 
1971. Using a realistic flow model which accounts for 
thermal and chemical nonequilibrium, he numerically obtained 
the flow characteristics for argon cascade arcs. The agree- 
ment of his rather extensive set of solutions with experi- 
mental data indicated that the flow model satisfactorily 
predicted arc flow behavior over a broad range of operating 
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conditions. For this reason, Clark's flow model has been 
selected for use in this study. 

Spectroscopic studies of hydrogen cascade arcs have 
been done by various researchers and a brief summary of 
their work is appropriate, Morris and others [10] report 
temperature profiles and degree of nonequilibrium (percent 
difference in electron and gas temperature) for currents of 
20 , 30, 40, 50 and 70 amps at one atmosphere pressure. They 
also report the electric field-current (E-I) relationship 
for a cascade arc of radius 0.0015 meters. Wiese and others 
[27] also report temperatures for a 40 amp atmospheric arc 
with a tube radius of 0.0015 meters. Measurements taken by 
Maecker [11] of the E-I relationship of a hydrogen arc with 
a tube radius of 0.001 meters are shown in his work on 
transport properties in high power arcs. Temperature pro- 
files and the degree of nonequilibrium for currents of 8,, 20, 
40, 60, 90, 120 and ISO amps in a 0.001 m. hydrogen arc have 
been obtained in the experimental studies of Steinberger [12] . 
Many of the above experimental results, as well as the 
numerical LTE studies of Greene, will be compared with the 
calculations of this study in a later chapter. 

1.4 Objectives of this Study 
It is the objective of this study to obtain solutions 
to a rigorous flow model of the hydrogen cascade arc that is 
not limited by unrealistic assumptions. As discussed in 
Section 2.1 on assumptions used in the flow model, only two 
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modes of nonequilibrium are found to be significant for a 
rigorous treatment of the hydrogen arc. These are the second 
and fourth modes of nonequilibrium which are mentioned in 
Section 1.2: nonequipartition of kinetic energy between the 

electronic and heavy particle species and chemical nonequi- 
librium of the various species present. The term ’’thermo- 
chemical nonequilibrium" will hereafter be used to describe 
this particular state of nonequilibrium in the arc. 

The primary purpose of this study is to obtain added 
insight into the general physics of the cascade arc by 
parametric studies with the hydrogen arc flow model, compari- 
son of the rigorous nonequilibrium solutions with equilibrium 
and experimental results for hydrogen arcs, and by comparison 
of the nonequilibrium hydrogen and argon [9] numerical data. 

This work is a continuation of an extensive theoretical 
and experimental study of nonequilibrium cascade arc flow at 
the Purdue High Temperature Gas Dynamics Laboratory. 
Theoretical work on thermochemical nonequilibrium in argon 
has been completed by Clark. This study extends the knowl- 
edge of nonequilibrium arc phenomena and the theoretical 
methods of treating it. 
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CHAPTER 2 MATHEMATICAL FLOW MODEL 
2 . 1 Assumptions 

In this section the assumptions made in the derivation 
of the governing equations are listed and discussed and in 
Section 2.2 the derivation itself appears. 

The assumptions which follow are numbered for future 
reference . 

1. The flow field is steady and all properties are sym- 
metric about the axis (axisymmetric) . At the flow 
rates considered in this problem, the flow field is 
laminar. 

2. Gravitational effects are negligible. 

3. Externally applied magnetic fields are absent. For 
argon, Bower [5] found that induced magnetic fields 
were negligible for currents less than 1000 amps. 

That assumption is assumed to be equally valid for 
hydrogen. 

4. The derivation of the macroscopic governing equations 
is valid from the continuum point of view since all the 
particle mean free paths are several orders of magni- 
tude smaller than the tube diameter. 
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5. Gradients of flow properties in the axial direction 
are much smaller (and, in fact, for the asymptotic 
region are identically zero) than those in the radial 
direction. Therefore, the usual boundary layer 
assumptions are valid, i.e. axial diffusion of heat, 
mass and momentum are negligible and the pressure is 
uniform over the arc cross section. This assumption 
is further justified by Clark [9] who utilized it in 
his nonequilibrium studies and obtained good agreement 
with experimental results. 

6. The governing equations are applied in regions suffici- 
ently far removed from the cathode and anode that the 
electric field has a nonzero component only in the 
axial direction. From Maxwell's equations this 
further implies that the electric field is uniform 
over the arc cross section. 

7. Flow rates corresponding to Mach numbers significantly 
less than unity are anticipated; therefore shear and 
pressure work and fluid kinetic energy are negligible. 

8. In the hydrogen arc, there are four chemical species 

that are present in significant amounts under equi- 
librium conditions. They are: electrons, positively 

charged atomic ions, neutral atoms and neutral molecules. 
At equilibrium, the concentration of all other species 

is less, by several orders of magnitude, than those 
just mentioned and it is assumed that these concentrations 
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remain insignificant in the nonequilibrium case as well. 
In this study, the presence of the molecular species 
has also been ignored for reasons which are discussed 
in Chapter 4. Therefore, governing equations are 
derived considering the presence of electrons, positive 
ions and neutral atoms. 

9, Strictly speaking, Maxwellian velocity distributions 

do not exist in a gas wherever gradients in flow proper- 
ties exist. However, for most problems the deviation 
from a Maxwellian distribution is slight enough to be 
negligible (and so the terms Maxwellian and near- 
Maxwellian are used interchangeably) . As mentioned 
in Chapter 1, Maxwellian velocity distributions are 
assumed for each of the chemical species. Incropera 
and Viegas [1] have found that in an argon cascade arc, 
this assumption is valid except near the cool wall. 
However, their results have been obtained using an 
equilibrium flow solution which significantly under- 
predicts the electron temperature and number density 
near the wall; hence the electron self-collision fre- 
quency (which is the principal equilibrium restoration 
mechanism) is substantially underestimated. Therefore, 
it is reasonable to expect Maxwellian velocity distri- 
butions to exist in the argon cascade arc. Since the 
mechanisms which promote nonequilibrium velocity dis- 
tributions in argon are the same as those in hydrogen, 
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Maxwellian distributions are assumed to exist in 
hydrogen as well. 

10. Equipartition o£ translational energy in a multicom- 
ponent mixture of gases occurs by means of elastic 
collisions between the particles of the different 
species. It is well known that, if two particles 
of approximately equal mass collide elastically, the 
exchange of energy between those particles is highly 
efficient. However, if the two particles differ 
greatly in mass, then the energy exchange process is 
highly inefficient and each particle leaves the colli- 
sion with practically the same energy it had prior to 
collision. Therefore, when two groups of particles 
of nearly equal mass but different thermal speeds are 
mixed, their translational energy equilibrates rapidly. 
This is not the case, however, with particles of greatly 
different mass. As far as this problem is concerned, 
it is therefore reasonable to assume that all the 
species consisting of heavy particles (ions and atoms) 
share a common mean energy and temperature, but since 
the electron is more than three orders of magnitude 
lighter than a heavy particle and is selectively 
energized by the electric field, the electron kinetic 
temperature cannot be assumed equal to the heavy 
particle temperature. Therefore, in the formulation 
of the flow model, it is assumed that all heavy particle 
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species have a common temperature and that the electron 
species has a different temperature. 

11. The perfect gas equation of state is assumed to apply 
to each chemical species and Dalton's law of partial 
pressures is assumed valid. Griem [25] calculates a 
correction for the perfect gas equation in a plasma 
due to Coulomb interactions and concludes that it is 
negligible for conditions of this study. 

12. In their discussion of electrical neutrality, Holt and 
Haskell [13] have shown that for ionized gases at the 
temperatures and electron number densities encountered 
in this work, conditions are electrically neutral on a 
macroscopic scale. For this work, the implication is 
that the hydrogen gas is electrically neutral every- 
where except in a microscopic region termed the plasma 
sheath (discussed in Section 2.2.3) which adjoins the 
wall. 

13. Due to large gradients in the electron and ion number 
densities, these particles will tend to diffuse from 
the core of the arc toward the wall. The condition 
of charge neutrality implies that, in the absence of 
externally applied forces in the radial direction (such 
as in this case) , the electrons and ions must diffuse 
together toward the wall. This condition of equal 
electron and ion diffusion velocities, or ambipolar 
diffusion, is considered in the derivation of the 
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governing equations . This assumption breaks down in 
the electrostatic sheath region which is discussed 
in Section 2.2,3. 

14. As a result of the large value of the electronic energy 
for the first excited level of atomic hydrogen, the 
population of this level will always be quite small, 
even in the core of the arc. A rough calculation 
(based on equilibrium population densities) shows that 
the maximum electronic internal energy contribution 
from the first excited level in the core of the arc is 
less than 1% of the particle kinetic energy. Therefore, 
electronic excitation effects are completely ignored in 
the consideration of energy storage mechanisms. Chemi- 
cal equilibrium, however, is not assumed to exist any- 
where in the arc, and, in fact, the treatment of the 
complications arising from this form of nonequilibrium 
is an important feature of this work. 

15. The presence of a radiation field, if treated rigorously, 
increases the complexity of the governing equations. 

For this reason, the hydrogen gas is assumed to be 
optically thin, and all radiation emitted by the gas is 
assumed to reach the cascade arc wall uninhibited by 
absorption. This assumption is not entirely justified, 
of course, due to the fact that the resonance radiation 
(that due to spontaneous transitions to the ground 
level) is quickly absorbed as it propagates through 
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the arc. Proper treatment o£ the radiation problem 
would involve the inclusion of the radiative transfer 
equation in the set of governing equations, a task 
which is well beyond the scope of this work. 

16. In order to calculate certain terms that appear in the 
governing equations, it is necessary to consider the 
microscopic nature of the gas. Transport properties, 
ionization rates and other properties depend on the 
collision frequencies of the various particles among 
and between each other. Therefore, the collision fre- 
quency is an important piece of information. However, 
in some high temperature gases, the notion of a colli- 
sion is vague since the trajectory of a particle might 
be simultaneously influenced by more than one neighbor- 
ing particle. However, Delcroix [14] points out that 
for conditions of interest here, this is not the case 
and most collisions can be treated as binary. The 
exceptions are those collisions involving the recom- 
bination of ions and electrons in the presence of a 
third particle (which receives the recombination 
energy) , the inverse of which is a two-body (ioniza- 
tion) collision. 

The potentially most important collisions which 
need be considered in this analysis are listed below. 
First the elastic collisions, for which the colliding 
particles suffer no change in kinetic energy, are 



22 


listed and then the inelastic collisions, which involve 
a conversion of kinetic energy to some other form, are 
given . 
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All of the collisions listed under Elastic Colli- 
sions are considered in deriving the transport proper- 
ties of hydrogen. All of the inelastic collisions were 
considered in the derivation of the source terms of the 
governing equations except atom ionization because data 
on the atom ionization rate could not be found. A 
study of atom-atom excitation in hydrogen [29]''indicates 
that atom ionization may be important in the wall region 


t 
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of the arc. However, the uncertainty analysis in 
Chapter 4 indicates that when the electron ionization 
(and recombination) coefficient is perturbed by a factor 
of 2, the effect on the solution is completely negligible. 
Therefore, it is expected that negligence of atom ioniza- 
tion has no effect on the solutions even if the atom 
ionization rate is as large as that for electron 
ionization. 

2 . 2 Derivation and Discussion of Governing Equations 

2.2.1 The Governing Equations 
Employing the assumptions discussed in Section 2.1, the 
governing equations are developed for the flow of an atomic 
gas through the entrance region of a cascade arc. The equa- 
tions for the asymptotic region are then a special case of the 
entrance region equations. Detailed derivations are not 
given, however, the physical significance of each term in 
the equations is indicated. 

In order to derive the necessary equations, the laws of 
conservation of mass and energy and Newton’s Second Law must 
be applied to the gas flowing through a cylindrical differ- 
ential volume fixed inside the cascade arc. A diagram of a 
cylindrical control volume is shown in Figure 2-1. 

The general form of the species continuity equation is 
derived by applying the conservation of mass principle to 
some species i as it flows through the control volume. 
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Volume 

Figure 2-1. Infinitesimal Cylindrical Control Volume. 
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Utilizing the appropriate assumptions from Section 2.1, this 
equation then becomes a balance between the net mass of i 
convected out of the control volume and the rate of creation 
of mass of i due to chemical reactions. 

Species Continuity Equation 

(1) (2) C3) 

net rate of change net rate of change net rate of change 

of species i in of species i in of species i in 

control volume due control volume due “ control volume due 

to axial convection to radial convec- to chemical reaction 

tion 

(1) (2) (3) 

3[p. (U+U.)] ^ ^ 3[PiCV+V^)r] _ . 

__ + _ _ 

Due to the assumed lack of large temperature and concen- 
tration gradients in the axial direction, the only diffusion 
force acting in that direction is due to the electric field. 
Therefore for neutral atoms is zero. Furthermore, due to 
the large relative mass of ions, the diffusion velocity of 
ions is small compared to that for electrons. Therefore, 
the only significant axial diffusion velocity is that of the 
electrons. This velocity will be referred to as the electron 

drift velocity . Discussion of radial diffusion veloci- 

* 

ties and the source terms Cpj^) is presented in Section 2.2.2. 

Applying the conservation of mass principle to the over- 
all gas mixture flowing through the control volume gives the 
following. 
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Global Continuity Equation 


( 1 ) 


( 2 ) 


(3) 


net rate of change 
of mass in control 
volume due to 
axial convection 


net rate of change 
^ of mass in control 
volume due to 
radial convection 


net rate of change 
of mass in control 
volume due to 
chemical reactions 


( 1 ) ( 2 ) (33 




^ 1 8[pVr] 


0 


Since total mass cannot be created or destroyed, term (3) is 
zero . 


Since the usual boundary layer assumptions have been 
made, the global radial momentum equation is not needed. 

The global axial momentum equation is derived by applying 
Newton's Second Law to the fluid flowing through the control 
volume . 


Global Axial Momentum Equation 


(13 


( 2 ) 


(3) 


resultant of all 
axially directed 
forces acting on = 
fluid in control 
volume 


net rate of change 
of axial momentum 
in control volume 
due to axial 
convection 


net rate of change 
of axial momentum 
+ in control volume 
due to radial 
convection 


(1) (2) (3) 


, dP 1 

dz" r 9r J 


9[pUM + 1 S[pUVrl 
3z r 3r 


The method of determination of the shear stress is discussed 
in Appendix B. 



27 


The electron energy equation is derived by applying the 
law of conservation of energy to the electron subgas flowing 
through the control volume. 


Electron Energy Equation 

(1) (2) (3) (4) 

net rate of net rate of net rate of net rate of 

change of elec- change of elec- change of elec- change of elec- 
tron energy in tron energy in _ tron energy in ^ tron energy in 

control volume control volume control volume control volume 

due to axial due to radial due to heat due to volumet- 

convection convection conduction ric sources 

and sinks 


( 1 ) 


( 2 ) 


5z r 3r 


(3) 

l3("Sre] 


(4) 


3r 


e ec icj 


In term (4), the electrical heating term is a 

source since it always imparts energy to the electron subgas. 
The elastic collision term (S ) is always a sink when the 
electron temperature is higher than the heavy particle tem- 
perature. The inelastic term (S®^) is neither source nor 
sink throughout the entire arc region but is a sink where 
ionization is predominant and a source when recombination is 
predominant. The term S? includes the radiation emitted 

At ^ 

by the gas. 

The electrical conductivity expression is derived in 

Appendix C and the electron heat conduction term (q ) is 

^e 

obtained in Appendix B. The remaining source terms (S and 
sf ) and the diffusion velocities (V and U,) are discussed 

1C- 6 CL 
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in Section 2.2.2. Derivation of the diffusion coefficients 
is given in Appendix A. 

The heavy particle energy equation is derived by apply- 
ing the law of conservation of energy to the heavy particle 
species (ions and neutral atoms) flowing through the control 
volume . 


Heavy Particle Energy Equation 


( 1 ) ( 2 ) 

net rate of net rate of 

change of change of 

heavy part- heavy part- 

icle energy ^ ical energy 
in control in control 

volume due volume due 

to axial to radial 

convection convection 


(3) (4) 

net rate of net rate of 

change of change of heavy 

heavy part- particle energy 

ical energy in control vol- 

in control ume due to vol- 

volume due umetric sources 

to heat and sinks 

conduction 


( 1 ) ( 2 ) 
3lPiUh.]j^ , 3[Pir(V.V.)h,]j^ 


dz 


3r 


(3) 

1 

r Sr 


(4) 


t- See) 


The bracketed terms [ indicate a summation over the 

heavy particle species. For example, ° 

The source term (- S ) is of the opposite sign as that in 

the electron energy equation since the heavy particles are 

the recipients of the energy lost by the electrons in elastic 

collisions. The term q_ is the contribution of the i^^ 

i 

species to the radial heat conduction. Derivation of the 
expressions for q_ is given in Appendix B. 

The equation of state concludes the list of formal 
governing equations that mathematically describe this flow 
problem. 
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Equation of State 


P 




• { OlVh' 


At the conclusion of this subsection the governing 

equations for flow in both the entrance and asymptotic 

regions of the cascade arc are listed for easy reference. 

In order to obtain the set of equations modeling the flow 

in the asymptotic region from the entrance region equations, 

the following steps are taken. The mass average radial 

velocity (V) is set equal to zero and all axial derivatives, 
dP 

except are also set equal to zero. Looking at the 
entrance region equations it is seen that this process 
eliminates the global continuity equation by forcing both 
terms to be equal to zero and the momentum equation is no 
longer needed since the axial velocity (U) does not appear 
in any of the remaining equations. The asymptotic set of 
governing equations then contains two fewer equations than 
the set of equations governing the entrance region, 
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Governing Equations of the Entrance Region 
Electron Continuity Equation 


3[p^(U+U^)] j 

5z r 


3[Pe(V+Vg)r] 
F 


* P, 


C2-1) 


Global Continuity Equation 



( 2 - 2 ) 


Axial Momentum Equation 


3Lp_uJJ 

3z 


1 3[pUVr] 
r 3r 


dP ^ 1 

Zz F “Tr 


(2-3) 


Electron Energy Equation 


§z r 5r r 


9[rq, 3 


3r 


(2-4) 

e ec 1 C 


Heavy Particle Energy Equation 


3[PlUhilh ^ 1 3(Pir(V*Vi)h.]^ _ 


3z 


3r 


3[rq_ 3 
1 ^i 
r 3r 
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Equation of State 
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(2-6) 
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Governing Equations of the Asymptotic Region 
Electron Continuity Equation 
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F §r 


(2-7) 


Electron Energy Equation 
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Heavy Particle Energy Equation 
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Equation of State 
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2.2.2 Properties Used in Governing Equations 
The dependent variables appearing in the entrance 
region governing equations are U, V, , h^ , P and E C^xial 
velocity, radial velocity, mass density of i, specific 
enthalpy of i , pressure and electric field) . It should be 
remembered that the above list of variables does not form a 
linearly independent set since the ionic and atomic enthal- 
pies are essentially equal and the electron and ionic mass 
densities differ by the electron- ion mass ratio. Hence, 
the dependent variables for this problem can be reduced to 
the following set; U, V, n^ , n^, T^, Tj^, P and E, where 

and n„ are the electron and atom number densities and T and 
H e 

are the electron and heavy particle temperatures. 

Each of the inelastic collision processes (listed in 
Section 2.1) can be quantified by the definition of a reac- 
tion rate, such that the product of the number densities of 
the reactants with the reaction rate coefficient gives the 
number of such reactions occurring per unit volume and time. 
For example, the rate of creation of electronic mass per unit 
volume due to the electron ionization of atoms could be given 
by: 

Pe = "‘e^e'^H^eH 

Therefore, in order to evaluate the term p appearing in the 
electron continuity equation, it is necessary to write it as 
the sum of those rate coefficients, multiplied by the 
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appropriate particle mass and number densities, which involve 
a reaction resulting in the production or consumption o£ an 
electron. This therefore requires obtaining the rate coeffi- 
cients for the forward and reverse reactions of inelastic 
collision number one (Section 2.1) and of the forward reac- 
tion of collision three. (Recall that collision two is 
insignificant.) For this study, the coefficient for colli- 
sion three was obtained from Allen [15] . The coefficient 
for the reverse reaction of collision one (electron- ion 
recombination) was obtained by fitting the data of Hinnov 
and Hirschberg [17] with a simple analytical function. The 
ionization rate coefficient is then obtained by using the 
law of detailed balancing together with the recombination 
rate coefficient. 

A 

The term, appearing in the electron energy equa- 

tion, is closely related to the source term in the electron 
continuity equation. Whenever an electron ionization (or 
recombination) collision occurs, there is a loss (or gain) 
of energy in the electron subgas in the amount of the ioniza- 
tion potential. Therefore, the product of the appropriate 
number densities with the rate for the forward (or reverse) 
reaction of inelastic collision one (Section 2.1) and the 
ionization potential of hydrogen gives the rate of loss (or 
gain) of energy by the electrons due to the reaction. As 
for radiative recombination (inelastic collision three) , 
there is a loss of energy in the amount of the ionization 
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potential plus the mean electron energy (1.5 kTe) when such 
a reaction occurs, and the contribution to is easily 
obtained. The Remaining contribution to enters as a 
loss of electron kinetic energy due to Bremsstrahlung 
(inelastic collisions four and five). The derivation for 
the energy losiJ rate due to electron- ion Bremsstrahlung is 
given by Clark [9] and that for electron-atom Bremsstrahlung 
is found in Appendix E. 

The electron drift velocity (U^) appearing in the 
electron continuity and energy equations is defined as 
follows [13] 



where e is the absolute value of the electronic charge and 

o is the electrical conductivity. The derivation of the 
© 

electrical conductivity expression is given in Appendix C. 

The transport properties (diffusion coefficients, 
viscosity, thermal conductivity) needed in this study, are 
generally provided in the literature in terms of their 
equilibrium values. Since the most important aspect of this 
work is that it treats nonequilibrium flow, equilibrium 
transport coefficients cannot be used. Therefore, a rela- 
tively accurate application of mean-free-path theory has 
been used to obtain these quantities. The derivations of 
the expressions for the ambipolar diffusion velocity (V^) 
and coefficients, the shear stress (t ) and the electron 

* lU 
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and heavy particle heat fluxes (‘Ij.g found in 

Appendices A and B. 

Finally, in a report by Petschek and Byron [18], an 
expression is derived for the rate of energy loss by elec- 
trons due to elastic collisions with a heavy particle subgas. 
Their expression is used in this work in order to obtain the 
equation for 

2,2.3 Boundary Conditions 

Since the boundary conditions of the asymptotic equa- 
tions are a subset of those of the entrance region equations , 
the boundary conditions of the entrance region are developed 
in this section. Specific discussion of the asymptotic 
boundary conditions is given where appropriate. 

In the governing equations of the entrance region, the 
axial velocity (U) , mass densities (p^) and specific enthal- 
pies (h^) appear in second order radial partial derivatives 
when appropriate substitutions are made for the shear stress, 
ambipolar diffusion velocity and heat flux, respectively. 
Therefore, two boundary conditions must be supplied for each 
of these variables. For this problem, conditions are speci- 
fied at the tube centerline and wall. Since the radial 
velocity appears only in a radial first derivative, its value 
needs to be specified at only one boundary (in this case, the 
centerline). Note that since E and P are not differentiated 
with respect to r, no boundary conditions are required. The 
centerline conditions will be discussed first. 
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Due to the assumed axial symmetry of the flow field, 
all radial transport phenomena must vanish at the tube 
centerline. Therefore, V Cranial velocity) and the radial 
first derivatives of U, and h^ are all zero at the tube 
centerline. The finite-difference version of the derivative 
boundary conditions yields the result that the values of U, 
and h^ at the centerline equal those at the first point 
away from the centerline (Figure 3-1) . However, a more 
rigorous set of centerline conditions may be obtained by 
utilizing the governing equations in conjunction with the 
first derivative boundary conditions. The governing equa- 
tions must be first evaluated in the limit (using L'Hospital’s 
Rule) as r and V approach zero, after which the boundary con- 
dition may be substituted. The "centerline” governing equa- 
tions, derived according to the above procedure appear on 
the following page. Note that all radial first derivatives 
appearing in the ambipolar diffusion velocity, shear stress 
and heat flux terms are zero in the centerline equations. 

Since the value of V has been specified at the center- 
line, it remains to specify boundary conditions on U, p^ 
and h^ at the wall. The usual nonslip flow condition is 
imposed at the wall, and hence, the axial velocity at the 
wall is zero. Since the atom number density is between three 
and four orders of magnitude larger than the electron number 
density at the wall and the atom-atom mean free path is 
roughly 5x 10"® centimeters, it is reasonable to assume that 
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Centerline Boundary Conditions of the Entrance Region 


Electron Continuity Equation 


3[p (U+U,)] 3(V+V ] . 

= p, 
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Heavy Particle Energy Equation 
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the heavy particle subgas will equilibrate with the wall. 
Therefore, Tj^ at the wall is equal to the wall temperature, 
which is fixed at 1000®K. This is equivalent to specifying 
hj_j+ and h^ at the wall. The atom mass (or number) density 
at the wall is determined from the equation of state and 
knowledge of the wall conditions for the remaining densities 
and enthalpies. It now remains to establish conditions for 
the electron temperature and number density at the wall. 

Due to the nonequilibrium assumption, the electron tem- 
perature and number density at the wall are unknown. How- 
ever, the necessary conditions may be obtained from considera 
tion of what is known as the wall sheath phenomenon. Since 
the electron thermal speed is much higher than that of the 
other heavy particles, the collision frequency of electrons 
with the wall is higher than that of the ions. This causes 
the metallic wall (which does not carry a current) to main- 
tain a negative charge relative to the gas in the constrictor 
The wall sheath region extends from the wall to a point where 
the gas becomes electrically neutral (about ten Debye lengths 
away from the wall) . By deriving conservation equations for 
the flux of electrons and electron energy across the sheath, 
the electron continuity and energy equations can be employed 
at the wall to yield two more boundary conditions. 

Clark [9] gives a detailed derivation of the wall heat 
flux and ambipolar diffusion velocity at the edge of the 
sheath region, therefore, only the expressions for these 
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quantities will be given here. The interested reader is 

referred to the work o£ Clark [9] and Knight [16] . 

The outward radial flux of electrons at any point in 

the gas is given by the product of the electron number 

density and the ambipolar diffusion velocity. The ambipolar 

diffusion velocity at the wall does not vanish because of 

electron-ion recombination at the wall. If 3 is the fraction 

(0 <. 3 <. 1) of electrons reflected from the wall without 

recombination, the flux of electrons to the wall is given by 

the product of the wall electron collision frequency, fg , 

w 

and the quantity (1 - 3). Therefore, the following equality 
holds 

ne V = f (1 - 3) C2-17) 

w w w 


where both sides of the equation express the flux of elec- 
trons to the wall. Equation (2-17) is the expression used 
for the electron flux in prescribing the wall boundary condi- 
tion for the electron continuity equation. Clark [9] derives 

the expressions needed to write f as a function of wall 

®w 

dependent variables . 

An energy balance for the sheath region yields the 
expression 


(1) (2) (3) (4) (5) 

(v K * I kVf = 2kT e|*J 

e w w w w 


(2-18) 
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where the terms in the equation are 

1. Fourier heat conduction into sheath from gas, 

2. Ambipolar diffusion of electron energy into sheath 
from gas, 

3. Energy transfer to the sheath due to electron reflection 

from the wall , where it is assumed that the reflection 
of electrons is diffuse and that the reflected elec- 
trons possess a thermal energy distribution character- 
istic of the wall temperature [16] . Therefore the 
transfer of energy from the wall into the sheath is 
given by the product of | with 6 and the collision 

frequency of electrons with the wall. 

4. Clark [9] shows that the average electron arriving at 

the wall has an energy equivalent to 2kT . Therefore, 

w 

the rate of energy transfer from the sheath to the wall 

is given by the product of the wall electron collision 

frequency and 2kT as in term (4) . 

®w 

5. Finally, every electron which successfully traverses 
the sheath potential to reach the wall and recombine 
with an ion loses energy in the amount e|<ti^|, where 
is the electrostatic wall potential. Term (5) gives 
the total loss due to this phenomenon. 
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Note that in Equations (2-17) and (2-18) several o£ the 
terms should be evaluated at the sheath edge rather than the 
wall. However, since the sheath is about ten Debye lengths in 
thickness, the difference is indistinguishable. Equation 
(Z-IS). When solved for then provides the wall condl- 

tion needed to solve the electron energy equation. The 
preceding wall boundary conditions appear on the following 
page for ease of reference. 

The centerline and wall boundary conditions of the 
asymptotic region are obtained from those of the entrance 
region by simply ignoring the conditions on the axial and 
radial velocities. The axial momentum and global continuity 
equations are no longer necessary, and the remaining equations 
may be simplified by setting all radial velocities and axial 
derivatives equal to zero. It should be noted that, in the 
computer program used to obtain solutions in the entrance 
region, the sheath conditions were deleted for simplicity, 
and conditions were obtained by simply specifying the elec- 
tron temperature and number density at the wall. 

No further conditions are required for the asymptotic 
equations, but since there appears in the entrance region 
equations, first order axial derivatives of U, and h^, 
these quantities must be specified at the constrictor 
entrance (z = 0) and are referred to as ’’initial conditions." 
The initial condition on the axial velocity has been chosen 
to be a parabola with a centerline value that satisfies mass 
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Wall Boundary Conditions of the Entrance Region 
Axial Velocity 

U(R) = 0 

Heavy Particle Temperature 

Th(R) - ^wall ■ “'>0°'' 

Atomic Mass (or Number) Density 

P - PeRefe * I 

P = n^kT^ ♦ (n^ ♦ n„)kT^ 

Electron Wall Flux 
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flux requirements. Both temperature profiles are also 
assumed to be parabolic, with preselected centerline and 
wall temperatures. The electron number density profile is 
determined by specifying a third temperature profile and 
using it to compute the electron number densities from the 
Saha equation. The centerline and wall temperatures of 
this third temperature profile are adjusted to yield reason- 
able centerline and wall electron number densities. Finally, 
the entrance atom number density profile is evaluated from 
the equation of state and the other initial conditions. The 
assumed initial conditions are summarized on the following 
page. Note that CL stands for centerline, w for wall and 
R is the tube radius . 
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Initial Conditions of the Entrance Region 
Electron Temperature 


T = T + 


T - T 


®CL % 

\ ^ 

■ ST 

1 


Heavy Particle Temperature 


Ti, = T, + 


Tv, - \ 

^CL 


1 - 


m 


Electron Number Density 
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1 -IpAT’ 
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f2irm^kT» 
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P e 
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Atom Number Density 
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Axial Velocity 
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CL 


1 - 
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^CL ~ ^ 
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1 ^ rdT 
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ra = mass flux. 
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CHAPTER 3 NUMERICAL SOLUTION TECHNIQUES 
3 . 1 General Discussion 

The primary purpose o£ this study is to obtain accurate 
information concerning the flow of hydrogen through a cascade 
arc by using theoretical procedures. To obtain accurate 
results, it is necessary to incorporate numerous microscopic 
and macroscopic phenomena in the governing equations, thus 
causing the equations to assume a complex form. Expressions 
for the source terms and the diffusion velocities are compli- 
cated and involve three or four dependent variables. Due to 
the high degree of coupling and the presence of highly non- 
linear terms in the governing equations, exact solutions to 
the equations cannot be obtained. The only recourse is then 
to use finite-difference solution techniques. 

Finite-difference methods in general have several 
properties in common. The most basic of these common ele- 
ments is the idea of breaking up the real, continuous solu- 
tion domain into a finite set of discrete points. The solu- 
tion of a problem by finite differences therefore consists 
of calculating values of the dependent variables at these 
points. This is accomplished by replacing the derivatives 
in the governing equations by approximations involving the 
values of the dependent variables at neighboring points of 
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the finite-difference grid. Two important aspects of the 
procedure include determination of the spacing of the dis- 
crete solution points (grid, stencil or star) in the domain 
of the independent variables and the selection of the finite 
difference operators which will replace the derivatives. 

It is not the purpose of this chapter to present in 
great detail the particular finite-difference techniques 
utilized in this study but rather to present an overview, 
such that the reader can supply the details and derivations 
necessary for full understanding. Also, the numerical 
methods discussed in this chapter apply only to the solution 
of the asymptotic equations. The solutions to the equations 
of the entrance region require different techniques and are 
of relatively minor importance in this study. The reader 
who is interested in entrance equation solution techniques 
is referred to the work of Clark [9]. 

There are two basic choices available regarding the 
spacing of finite difference grid points: uniform or non- 

uniform spacing. Uniform spacing is to be preferred because 
it permits the simplification of the finite-difference 
operators and if, as in many cases, computer run time is 
small, the grid spacing can be easily reduced for greater 
accuracy. However, if computer time is expected to be large, 
as in this case, an optimally spaced variable grid is the 
wisest choice. The spacing can be made coarse in regions of 
small gradients and fine in regions of large gradients. 
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Therefore, the nonuniform grid, with accompanying compli- 
cated expressions for difference operators, has been selected 
for this work. The grid consists of 20 points with center- 
line spacing six times greater than at the wall. A schematic 
drawing of this grid appears in Figure 3-1. 

At first glance, it appears that all derivatives in the 
asymptotic equations (Section 2.2.1) are of the first order. 
However, the ambipolar diffusion velocity (V and V.) and 

the heat fluxes (q and q ) involve first order derivatives 

’^e ^i 

of the dependent variables. Therefore, each derivative term 
of the asymptotic equations, after substitution and expan- 
sion, takes the form, ^ » where (f> is a 

dependent variable and a is a single term or product of terms 
which are a function of (f> and possibly other dependent 
variables . 

The derivation of a finite- difference operator for the 

above derivative takes place in two steps. First, an expres- 

3 f 

sion is found to approximate evaluated at the grid point n 
in terms of the values of f at the half grid points (n - |) 
and (n + i) , Note that f(n + ^) = [f(n) + f(n + l)]/2. 

Half grid points are used for accuracy and for reasons that 
will be obvious in succeeding paragraphs. Deriving expres- 
sions for f(n + I) and f(n - 5 ) by expanding a Taylor series 
about the point n and manipulating the resulting two equa- 
tions in order to cancel the lowest order error term gives 
the equation 
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'9f' 

^ 2Ar. f 

2 ] 

Ar+ Ar.' 

[tf\ 

Ar+(Ar. + Ar+) n+^ 

Ar, + Ar+ 

t - T J 

Ar. ■ St7^ 


n 


ArJ Vi 


where Ar. and Ar+ are the differences in the independent 
radial variable, as shown in Figure 3-1. 

Substituting ra(<j ) , . . . )|y foi" f Equation (3-1); 
setting 


and 


(H) 

[57J 

II) 


- (1) 

^n+i 


n+| 

Ar+ 


A - A 
^n ^n- 

n-1 

Ar. 


evaluating ||± 


by means of Equation (3-1) by replacing 


n 


n+^ with n+1, n-| with n-1, Ar_/2 with Ar. and Ar^/2 with 
Ar+ as shown below 


l±] = At- * f 1_ t(Ar+ Ar. 


3rJ^ Ar+ (Ar. + Ar+) ^n + 1 [Ar. + Ar+J [Ar. ' Ar+ 


n 


Ar 


Ar.(Ar. + Ar+) *^n - 1 ’ 


the following expression for the desired derivative is 
obtained 




2(Ar.|. - Ar_) 
Ar?(Ar. + Ar+) 


Substitution o£ the above approximation in the governing 
equations yields the desired finite-difference approximation, 
the treatment of which is discussed in the following section. 
The truncation error of Equation (3-2) is complicated, but 
the lowest order terms are written below. 



The largest error (which is of the first order) is reduced 
by selecting a variable grid in which grid spacing increases 
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or decreases gradually. If this is done, the truncation 
error is of the second order. 

The centerline boundary equations, Equations (2-113 
through (2-16) , must also be approximated by use of second 
order finite-difference operators. However, the grid can 
be assumed uniform at the centerline and the second deriva- 
tives approximated by difference operators using values of 
variables at the centerline and the first point on either 
side of it (the variables at the latter two points being 
equal) . Since the standard expression for the second order 
difference operator (which has a second order error) for a 
uniform grid is well known, it will not be written here. 

A special difference operator must be developed for 
Equations (2-7) and (2-8) which are used in conjunction with 
Equations (2-17) and (2-18) for wall boundary conditions on 
electron wall particle and heat fluxes. Since the expres- 
sions for the heat flux and ambipolar diffusion velocity at 
the sheath edge do not contain derivatives, these equations 
are only of the first order in the radial derivative. 
Therefore, a finite difference expression must be obtained 
which approximates the first derivative of a quantity <|» at 
the wall. This is accomplished by using a Taylor series 
expansion at the wall to approximate the value of the vari- 
able <|) at the points w-| and w-| (see Figure 3-1) in terms 
of the value of ^nd its derivatives at the wall. Two 
equations are obtained from this procedure and are solved 
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simultaneously in order to eliminate the second derivative 
terms in each. When this is done, the following expression 
is obtained for the first radial difference at the wall 


where 


3r 


w 


R6 3-Rd) i+Rd) 
Wj W-5 Wj w-2 ^3 w 


2Ar 


w 


w. 


f 2Ar + Ar ) (Ar + Ar T 
'■ w w - 1 w w - 1 


2(2Ar^ * 

2(3Ar^ + 

^w, ~ Ar ,(2Ar, , + Ar, , J 
3 w^ w w - 1 


and Ar and Ar are defined in Figure 3-1. The truncation 
w w - 1 ® 

error of the above formula is of the second order as shown 
below 

Ar (2Ar + Ar ) 
w^ w w-i-^ 




When this difference formula is directly applied to the wall 
electron continuity and energy equations, the resulting 
difference expression will contain the ambipolar diffusion 
velocity and heat flux evaluated at w-|, w-^ and w. At the 
first two of these three points, the standard differential 
forms of the ambipolar velocity and the heat flux are sub- 
stituted (see Appendices A and B) but at the last point, w, 
the sheath expressions. Equations (2-17) and (2-18), are 


substituted. 
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3 . 2 Solution of Finite Difference Gove rning Equations 
in order to introduce the method used to solve the 
difference equations, an illustrative example is first 
described. Suppose that it is necessary to solve a system 
of two transcendental equations in two unknowns . Due to the 
nonlinearity and coupling of the various terms of the equa- 
tions it is impossible to obtain a direct solution and, 
therefore, an iterative technique must be utilized. Assume 
that the unknowns in the equations are labeled x and y and 
that the equations have been written such that all terms are 
on the left hand side. Therefore, the system of equations 
appears as 

f 1 (x,y) = 0 ; fz Cx,y) = 0 

The Newton-Raphson method is to be utilized to obtain the 
values of x and y which satisfy the above expressions. First, 
expand each of the above functions to the first order in a 
Taylor series about x and y. 


i i+ 1 i+ I 

r ■>' J 

= f 

1 1 
X ,y 


i 

i+i 

X - X 

+ iiL 


, 

1 

^ 1 

9x 


3y 


y 

J 


= f 

2 

i i 

+ ^^2 
3x 

i 

f i+i i' 

X - X 

L J 

+ 3f2 

3y 

[y - 

J 


The above equations indicate that, if values of fi^* and f^^^ 
are arbitrarily selected, it is possible to find the set of 
variables (x,y) which will satisfy the above equations. It 
is desired now to find the values of x and y which result 
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in fj and * being identically zero. Suppose a set o£ 
variables (x^,y^) have been arbitrarily substituted into the 
functions and have failed to yield zero. A second estimate 
of variables which will accomplish this is 

acquired by setting the left hand sides of the above equa- 
tions equal to zero and solving for Using 

matrix notation, this results in the following expression. 

" Mi. ' 

3y 

i Ml " 

ay 

Since the Taylor series expansion was terminated after the 
first order, the above matrix equation is only approximately 
correct. Therefore, the variables (x^*\ are only 

approximations of the true solution of the two transcendental 
equations. In order to obtain the true solution, the vari- 
ables (x^"^\ must be used in place of the values (x^,y^) 

to completely reevaluate the right hand side of the above 
matrix equation. Then, another matrix inversion, multiplica- 
tion and subtraction must be carried out to obtain a yet more 
refined set of variables (say . This process is 

continued until sufficient accuracy is obtained in the solu- 
tion. This method can, in general, be used for large systems 
of equations where it is possible to invert matrices quickly 
by the computer. 
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By utilizing the difference operators of the previous 
section to obtain the finite difference forms of the govern- 
ing equations and the boundary conditions, a large system of 
equations is generated in which the unknown variables are 
the densities and enthalpies at each of the individual grid 
points. Actually, the final form of the governing differ- 
ential equations used in this study has been simplified so 
that the variables are electron and atom number densities 
and electron and heavy particle temperatures. Therefore, 
there are four unknowns in the four governing equations of 
the asymptotic region (Section 2.2.1). 

For each grid point at which it is desired to obtain a 
solution, there are four unknowns (two temperatures and two 
number densities) and four finite-difference governing equa- 
tions. Suppose there are n such grid points. Therefore, 
there are 4n unknowns (n , n„ , T , T. , i = l,n) and 4n 

Cf e^ 

finite difference equations. The Newton- Raphs on method has 
been employed in this study to solve this system of equations 
and its development proceeds as follows. 

The most direct approach is to choose an initial estimate 
of the variables at each grid point; calculate the values of 
the functions and their derivatives which appear in the 
matrix to be inverted, as in Equation (3-3), (it is assumed 
that all terms in the difference equations are on the left 
hand side) ; invert the matrix; calculate the new set of 
variables; and then repeat the cycle until convergence is 
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obtained. In practice, this procedure did not yield con- 
vergence due to the strong coupling between equations and 
the following modification was employed. The electron con- 
tinuity equation and equation of state were "decoupled" from 
the two energy equations in the finite-difference scheme as 
follows. Initial estimates of variables at all grid points 
were first obtained, and all temperatures were assumed to 
be fixed. The matrix inversion routine then only involved 
the electron continuity equation and the equation of state. 
Therefore, only the number densities were recomputed until 
convergence was obtained. Then, the number densities were 
fixed and the matrix inversion routine involved only the 
energy equations. The temperatures were then recomputed 
until convergence was obtained. This procedure was repeated 
until successive cycles produced unchanging number density 
and temperature profiles, indicating total convergence. It 
was also necessary to introduce the following additional 
step in the actual computer program. Rather than specify 
the electric field (which appears explicitly in the electron 
energy equation) , the current was actually read into the 
computer program and the electric field was estimated. 
Therefore, after convergence was obtained on the temperature 
profiles (with fixed number densities) , the integration of 
the electrical conductivity was carried out to calculate the 
current corresponding to the estimated electric field. This 
was done according to the expression 
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I = E I 27ra^rdr (3-4) 

If the calculated current did not agree with that which was 
read into the program, an iteration procedure (Newton- 
Raphson) was initiated to obtain the correct electric field 
by estimating a new value of E, substituting it back into 
the energy equations (with fixed number densities) and obtain- 
ing the new converged set of temperatures and electrical con- 
ductivity. Integration was again repeated according to 
Equation (3-4) and the process continued until the correct 
electric field was obtained. 

One final comment needs to be made regarding the selec- 
tion of the first estimate of the unknown variables. It has 
been found that this program is sensitive to the accuracy of 
the input variables. For this reason another computer pro- 
gram (program MARCH) was developed to use an explicit, 
finite-difference scheme to solve the governing equations of 
the entrance region. As the solution of program MARCH pro- 
gresses in the axial direction of the cascade arc, it 
approaches the asymptotic region and yields number density 
and temperature profiles suitable for the initial values of 
the asymptotic solution method. However, most of the initial 
profiles used in this work were provided by solutions of the 
asymptotic region program for other operating parameters 
(current, pressure, tube radius). Therefore, due to the 
insignificance of program MARCH, the explicit finite-difference 
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scheme for the entrance region equations has not been dis- 
cussed here. The interested reader is referred to the work 
of Clark [9] . 
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CHAPTER 4 RESULTS AND DISCUSSION 
4 . 1 Preliminary Discussion 

Development of the finite-difference and computer pro- 
gramming techniques in order to obtain stable, convergent 
solutions to the governing equations has proved to be a 
formidable task. It is not possible to mention all of the 
intricate details which have become necessary to obtain solu- 
tions, however, some of the more important aspects of the 
methods should be discussed. 

As mentioned previously, a computer program (MARCH) was 
developed to solve the governing equations of the entrance 
region. Since this program proved to be stable only at low 
currents, it was of limited value and has not been discussed 
to any great extent. However, this program has been used to 
obtain a solution for I = SO amps, P = 1 atm. and R = 0,005 
meters, and the results were sufficiently accurate to serve 
as the starting data for the iterative solution (the Newton- 
Raphson method) of solving the asymptotic equations. In 
this study, three separate computer programs were used to 
obtain the solution to the asymptotic equations in final 
form. The first program (NRHEQD) provided solutions to the 
equations for the case of thermal equilibrium (T^ = Tj^) with 
fixed (rather than electrostatic sheath) boundary conditions. 
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This program was used to obtain approximate densities and 
temperature for new operating conditions because it was 
faster than the other two. When a set of profiles was 
obtained from program NRHEQD, they were read into a second 
program (NRHD) which obtained the solution to the governing 
equations without the thermal equilibrium restriction but, 
again, with fixed boundary conditions. This then provided 
a sufficiently accurate set of data for the more sensitive 
program (NRHDS) which accounted for the sheath electrostatic 
wall phenomena. The output from NRHDS was then the final 
solution. The computer time required to obtain each solution 
by this procedure on a CDC 6500 computer ranged from about 
10 minutes to an hour. 

In order to maintain numerical stability of the asymp- 
totic computer programs at higher currents, it became neces- 
sary to fix the electron and heavy particle temperatures at 
the grid point adjacent to the wall (w-1 in Figure 3-1). 
Reasonable values were selected on the basis of experience 
with the solutions, and solutions for the lower currents 
which did not require this restriction indicated that tem- 
peratures adjacent to the wall are relatively insensitive 
to values of the reflection coefficients. 

The criterion used to determine whether a solution had 
converged sufficiently was that successive iterates of the 
densities and temperatures at all grid points must not vary 
by more than 0.00051. Therefore, if the temperature at a 
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certain grid point is about 10,000°K, at the time o£ con- 
vergence it changes less than O.OS^K from one iteration to 
the next. In addition, a global energy balance was also 
used to check the convergence of each of the final solutions. 
The amount by which this balance was not satisfied averaged 
about 2 . 5 % for all parametric solutions. Undoubtedly, this 
figure could have been reduced by substantially decreasing 
the grid spacing, but this would have resulted in prohibitive 
computer run times. Also, the accuracy demonstrated by the 
convergence criterion and the energy balance is considered 
to be much greater than that of the basic data (such as 
recombination coefficients, transport properties, etc.) used 
in the program (see Section 4.4). Therefore, it was deemed 
of little benefit to obtain greater numerical accuracy. 

Solutions were obtained for two purposes. First, 
parametric predictions were obtained to gain more insight 
to the physical nature of high temperature, nonequilibrium 
arcs. In addition, solutions were obtained specifically 
for the purpose of comparison with experimental data and 
equilibrium solutions- The parametric predictions obtained 
covered the following range of operating conditions: 

1. P = 1 atm, R = 0.005 meters, 50 £ I £ 200 amps. 

2. I = 200 amps, R = O.OOS meters, 1 £ P £ 10 atm. 

3. I = 200 amps, P = 1 atm, 0.005 < R £ 0.05 meters. 
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Predictions were obtained for comparison with the 
experimental results of Steinberger [12] , Morris and Rudis 
[10], Wiese [27] and Maecker [11] and with the equilibrium 
predictions of Greene [6]. The range of operating conditions 
used in these comparisons is as follows: 


1. 

P = 1 

atm , 

R = 0.0015 meters, 20 < 

I < SO 

amps . 

2 . 

P = 1 

atm, 

R = 0.0010 meters, 20 < 

I < 40 

amps . 

3. 

P = 1 

atm. 

R = 0.005 meters, 50 < I 

< 200 

amps • 


Before discussing the solutions to the atomic hydrogen 
nonequilibrium model it would be beneficial to briefly con- 
sider the individual terras in the electron continuity and 
energy equations. 

Figure 4-1 illustrates the two terms of the electron 
continuity equation which exactly balance each other. It 
is evident that electron production in the core due to 
ionization is balanced by ambipolar diffusion of electrons 
away from that region. As the wall is approached, n 
reverses sign indicating that electrons diffusing into the 
wall region are being consumed by radiative and three-body 
recombination processes. 

Figure 4-2 is plotted in' such a way that positive 
values indicate energy transfer to the electron gas and 
negative values indicate energy transfer from the gas. As 
expected, the Ohmic heating term (OgE^) supplies energy to 
the electron gas over the entire arc cross section. Also 
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Figure 4-1. Terms in the Electron Continuity Equation. 
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Figure 4-2. Terms in the Electron Energy Equation. 
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radiative emission depletes the electron gas energy, 

especially in the core. The term (Pj) refers to the diffu- 
sion of ionization energy. Since electron-ion pairs are 
being produced in and diffused away from the core, Pj is 
negative there. However, as the wall is approached, the 
electron-ion pairs recombine transferring ionization energy 
to electron third particles and Pj is positive. The term 
^ic electron energy equation (2-8) is 

equivalent to * (Pj.g^j * core also loses thermal 

energy through ambipolar diffusion and electron heat con- 
duction as shown in Figure 4-2, but as the cooler wall 
region is approached the electron subgas begins to gain energy 
due to these processes. Finally, note that the exchange of 
elastic collisional energy between electrons and heavy parti- 
cles (S ) behaves unexpectedly. In fact, the heavy particles 
transfer energy to the electrons in the core region (for 
these operating conditions) indicating that the heavy particle 
temperature is slightly higher (by about .05%) than the elec- 
tron temperature. The reason for the magnitude of S being 
so large for such a small difference in temperature is because 
is proportional to the product of the electron density and 
the electron-heavy particle collision frequency, both of which 
are much larger in the core than elsewhere. From Figure 4-2 
it is evident that the high rate of depletion of electron 
energy in the core by electron heat conduction and electron 
ionization (as evidenced by the Pj term) tends to depress the 
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electron temperature. Note that from Equation (2-9), the 

integral - - [ rS -dr] gives the heavy particle heat con- 
[ r j 0 ec j 

duction and that from Figure 4-2 this indicates that there 
is a net heavy particle conduction of energy to the core 
from the periphery. This is impossible for simple Fourier 
conduction, but the mean free path theory of this model 
predicts a large diffusion- thermo effect (heat transfer 
induced by a concentration gradient) which is responsible 
for this effect (the heavy particle concentration increases 
rather than decreases with increasing r) . Nevertheless, at 
the wall, Fourier conduction takes over and results in a 
positive transfer of heavy particle energy to the wall. The 
consequence of this diffusion- thermo effect, coupled with the 
depletion of electron energy by conduction and ionization, is 
that the heavy particle temperature can be elevated slightly 
above the electron temperature because there is a source of 
energy for the heavy particle subgas in the core. Neverthe- 
less, the effect is small resulting in only a few degrees 
temperature difference. 

The comparative studies are discussed in detail in 
Section 4.2 and the parametric studies are discussed in 
Section 4.3. The chapter is concluded with an uncertainty 
analysis and a discussion of the molecular hydrogen arc in 
Sections 4.4 and 4.5, respectively. 
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4 . 2 Comparisons with Other Data 
In this section, comparison o£ the numerical solutions 
is made with both experimental results and equilibrium pre- 
dictions obtained by other researchers. Morris and others 
[10] have published the results of their study for an atmos- 
pheric hydrogen arc with a radius of 0.0015 m and sample 
temperature measurements are shown' in Figure 4-3. The curve 
shown is for a current of 40 amps and its relationship to 
the numerical results of this study is representative of the 
other currents for which data were obtained. Therefore, 
only the 40 amp curves are shown. The Morris profile is 
much more "peaked" than the results of this study, producing 
a higher centerline temperature. Wiese [27] has also obtained 
spectroscopic measurements of a hydrogen arc at the specified 
operating conditions and his results are also shown on Figure 
4-3. In contrast to the Morris results, the Wiese data agree 
closely with the results of this study. Actually, Wiese 
reports temperature profiles determined by two different 
methods. One set of data is obtained from spectroscopic 
measurements of the continuum radiation intensity and the 
other by measuring the intensity of the Balmer line of hydro- 
gen. According to Wiese [27] however, the continuum measure- 
ments are the less reliable of the two (due to the contribu- 
tion of several molecular processes to the continuum emission 
and the presence of weak, highly broadened impurity lines)'. 

For this reason only the temperature profile resulting from 
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Figure 4-3. Comparison of Theoretical and Experimental 

Temperature Distributions for a 40 amp, .0015 m 
Radius Hydrogen Arc. 
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Wiese's measurement of the Balmer line intensity is shown in 
Figure 4-3. His other temperature profile (not shown) is 
more in agreement with that of Morris, who obtained the 
temperature from a measurement of the continuum intensity. 

Figure 4-4 provides the comparison of experimental 
temperature profiles for a 0.001 ra radius atmospheric arc 
with the solutions of this study. In his experiments, 
Steinberger [12] also used the Balmer line intensity to 
obtain his temperatures. Note that at 20 amps, the agree- 
ment between theory and experiment is rather poor. Moreover, 
for this current, the predictions reveal the existence of 
thermal nonequilibrium over most of the arc cross section. 

Much time has been spent in thoroughly checking the solution 
procedures, computer logic, and programs of this study and 
there is no apparent reason why the low current solution 
should be any less valid than the high current results. It 
is to be expected that, at lower currents, nonequilibrium 
effects extend over a greater portion of the arc and that 
accurate intensity measurements are more difficult to obtain 
because of arc instabilities and asymmetry. The agreement 
between Steinberger 's data for 40 amps and the solution from 
this study is good. 

Greene [6] has compiled a set of numerical solutions 
for the mathematical model of a hydrogen arc in thermochemical 
equilibrium and representative results appear in Figure 4-5 
along with the solutions of this study. These solutions are 
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Figure 4-4. Comparison of Theoretical and Experimental Tem- 
perature Distributions for a .001 m Radius 
Hydrogen Arc. 
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Figure 4-5. Comparison of Temperature Profiles for Equi- 
librium and Nonequilibrium Hydrogen Cascade 
Arcs . 
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for an arc of radius .005 m and pressure of 1 atmosphere. 

Tlie comparisons of this figure reflect the effect of thermo- 
chemical nonequilibrium in the arc. 

At first glance, agreement seems to be poor, but it must 
be remembered that Greene's [6] solutions are for an arc in 
local thermodynamic equilibrium. This means that the number 
densities calculated near the wall are several orders of 
magnitude lower than for the nonequilibrium case. Therefore, 
his electrical conductivity is much smaller and the Ohmic 
heating term (a E^) in his energy equation is very small in 
the wall region. (See Figure 4-2). Note that, from Table 
4-1, although the square of his electric field intensity is 
as much as twice the nonequilibrium value, his electrical 
conductivity is orders of magnitude smaller. Therefore, the 
equilibrium model, as expected, greatly underpredicts the 
temperature in the wall region due to its failure to include 
Ohmic dissipation effects in this region. Furthermore, 
since the electric field is inversely proportional to the 
integral of the product of the radial coordinate and the 
electrical conductivity across the tube radius, the electric 
field of the equilibrium arc is higher due to the smaller 
electrical conductivities in the wall region. The elevated 
centerline temperatures for the equilibrium arc tend to com- 
pensate for this effect on the electric field but fail to do 
50 completely because of the small value of the radial 
coordinate and the fact that electrical conductivity is 
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comparatively insensitive to temperature at the higher arc 
temperatures (Figure C-1). Therefore, the higher electric 
field provides more Ohmic heating at the centerline in the 
equilibrium case and hence produces the higher centerline 
temperatures shown in Figure 4-5. 

As the current increases, the equilibrium and nonequi- 
librium electric fields converge (see Table 4-1). This 
effect is reflected by the results of Figure 4-5. As the 
current increases, the centerline temperatures converge, and 
good agreement between the theories exists over a larger 
portion of the arc cross section. The difference between 
the temperatures in the arc periphery is due to thermochemi- 
cal nonequilibrium effects which persist at the higher 
currents . 

In addition to determining the experimental temperature 
profiles, Morris [10] and Wiese [27] also measured the 
electric field intensity, and the results are presented in 
Table 4-1. Agreement between the experimental and nonequi- 
librium theoretical electric field intensities is poor at 
low currents but improves with increasing current. The best 
comparison is with the data of Wiese [27] which also provided 
the best agreement for the temperature profiles. 

The reason for the poor agreement between theory and 
experiment is thought to rest largely with the theoretical 
method because of convergence errors in the numerical pro- 
cedure at small tube radii. Consider the error in the 
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Table 4,1 Comparison of Experimental, Equilibrium and 
Nonequilibrium Electric Field Intensity Data 
and Predictions for Various Arc Operating 
Parameters . 


P = 1 atm, R = .001 m 


I ( amp s ) 

E (volts/meter) 

Maecker [11] 

. . 

Nonequilibrium 

20 

12000 

5572 

40 

9800 

5261 

50 

9000 

— 


P « 1 atm, R = .0015 ra 


I (amps) 

E (volts/meter) 

Morris [10] 

Wiese [27] 

Nonequilibrium 

20 

9300 

— 

4040 

30 

7400 

— 

3779 

40 

6800 

5600 

3672 

50 

6200 

1 

— 

3617 


P = 1 atm, R = .005 m 


I (amps) 

E (volts/meter) 

Equilibrium [6] 

Nonequilibrium 

50 

1824 

1343 

100 

1468 

1200 

150 

1375 

1191 

200 

1348 

1223 

300 

1361 

— 

400 

1398 

— 
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second order finite-difference formula given in Chapter 3 
and repeated here for reference. 


error 


CAr+ 


Ar_) 


3 ar3 


It is of interest to see how the error is affected by reduc- 
ing the tube radius and grid spacing by some factor "a” 

(less than unity) and assuming that the arc operating param- 
eters are adjusted such that the profiles of a and ^ maintain 
the same relative shape. Let the new values of tube radius 
and grid spacing be denoted by a primed (’) superscript such 
that r’ = ar and Ar’ = aAr' . The new error is then 


(Ar 


T 

+ 


Ar^) 


(r * a) 


n 
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”3 S(ar) 3 



Ar.) 
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Therefore if the number and relative spacing of grid points 
with respect to each other is maintained the same and the 
values of a and remain the same at each of the grid points, 
the error is proportional to To be specific, the profiles 
of temperature and number density for the solution at I = 40 
amps and R = .001 ra are roughly the same shape as profiles 
obtained for a radius of .005 m (which is used for parametric 
studies) and a current of 200 amps. Therefore, the finite 
difference convergence error is ^ = 5 times greater at 
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I = 40 amps and R = .001 m than at I = 200 amps and 
R = .005 m. In order to obtain the same accuracy at 
R = .001 m as at R = .005 m then it is necessary to increase 
the number o£ grid points from 20 to 100. Since computer 
time is roughly proportional to the square of the number 
of grid points for this type of problem, 25 times as much 
computer time would be required to obtain this accuracy. 

This would then involve computer run times of several hours 
per solution, which is indeed prohibitive. The above calcu- 
lation of the tube radius dependence of the convergence error 
is born out by the fact that at the tube radius of .001 m, 
energy balances of the final solutions showed discrepancies 
as high as 151, while at R = .005 m the discrepancies are 
less than II. 

It is unfortunate that experimental studies have not 
been performed for larger tube radii where the convergence 
error of this numerical method is tolerable. Although the 
solutions of this model are not accurate at low tube radii 
(R = .001 m) , good accuracy is expected at all the tube 
radii used in the parametric studies of Section 4.3. 

Comparison of the equilibrium and nonequilibrium pre- 
dictions of the electric field intensity is also shown in 
Table 4-1 as well as Figure 4-11. As mentioned previously, 
the difference in these predictions is due to the failure of 
the equilibrium model to account for thermochemical nonequi- 
librium effects . 
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Comparisons between the experimental and nonequilibrium 
total wall heat flux exhibit the same trends as the electric 
field comparisons, since the wall heat transfer is equal to 
the product of the electric field intensity and current in 
the asymptotic region. The equilibrium and nonequilibrium 
predictions of total wall heat flux are compared in Figure 
4-12 for both hydrogen and argon. Because of the larger 
value of the field intensity, the equilibrium solution pre- 
dicts a larger value of the heat flux than the nonequilibrium 
solution for the hydrogen arc. 

Finally, a comparison between equilibrium and nonequi- 
librium predictions of the radiation heat flux fraction is 
shown in Figure 4-13, along with the corresponding predic- 
tions for argon. Before discussing Figure 4-13, it is use- 
ful to discuss the radiation model used in this study. 

In the nonequilibrium model of this study, it is diffi- 
cult to rigorously treat the radiation power loss. Such a 
treatment would require solution of the conservation equa- 
tions for the atomic energy level populations and the equa- 
tion of radiative transfer simultaneously with the governing 
equations of this study. Since the inherent computational 
difficulties are prohibitive, an approximation to the radia- 
tive loss term is made. It is assumed that all electrons 
recombining radiatively with an ion would, after recombination 
to an upper level, radiatively decay to the ground state. 
Likewise it is assumed that any electron recombining with an 
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ion collisionally would continue to de-excite collisionally 
to the ground state. If such is the case, the energy radia- 
ted for each radiative recombination would equal the ioniza- 
tion potential plus the mean electron thermal energy 

3 

(Ip + j . Therefore, the radiative emission rate would 

be equal to the product of the radiative recombination rate 
and (Ip + j plus whatever energy is released due to 

Bremmstrahlung . 

The radiative recombination coefficient, though it can 
be accurately calculated for atomic hydrogen, is not strictly 
appropriate for use in the calculations of this study. In 
the calculation of this property from atomic theory, it is 
assumed that the presence of other modes (collisional) of 
recombination and excitation have no effect on the radiative 
recombination coefficient. This assumption is good for 
plasma afterglows and atmospheric radiation studies because 
of the very low pressures involved, but in an atmospheric 
pressure arc the collisional modes of recombination and 
excitation are dominant. This means that, when a radiative 
recombination takes place to an energy level other than the 
ground level, the probability of an electron, colliding with 
that excited atom and re-ionizing or de-exciting it is 
significant. This is especially true in hydrogen since the 
ionization energy for atoms in the upper energy levels is 
roughly equal to the mean electron kinetic energy. There- 
fore, the effect on the radiative recombination coefficient 
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of having a collisionally dominant gas is to reduce it by 
some factor less than unity. Comparing the radiative 
emission rate predicted by the above model, in the equilibrium 
limit, with the results of Yos [21], indicates that this mul- 
tiplicative factor is roughly 0.4 for hydrogen. Therefore, 
in this study, the low pressure radiative recombination 
coefficient given by Allen [15] is corrected for high pres- 
sure collisional interactions by the factor 0.4. 

Looking again to Figure 4-13 in conjunction with Figure 
4-12, the equilibrium radiative flux at 50 amps exceeds the 
nonequilibrium flux by 351 but at 200 amps the discrepancy is 
reduced to about 41. The reason for the rather large dis- 
crepancy at 50 amps is evident in Figure 4-5 where the equi- 
librium and nonequilibrium temperature profiles are plotted. 
The equilibrium core temperatures far exceed those at non- 
equilibrium. Since the radiative flux is a sensitive function 
of temperature, it is expected that at low currents the 
equilibrium model will overpredict radiative flux. Discus- 
sion of the equilibrium and nonequilibrium temperature com- 
parisons has been given previously in this section. Although 
for argon the equilibrium and nonequilibrium radiative flux 
fractions do not agree well, the nonequilibrium predictions 
agree very well with available experimental argon data. The 
difference in the argon and hydrogen radiative flux fractions 
is discussed in the following section. 
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4 . 3 Parametric Study 

There are three critical parameters which have been 
varied in order to examine the phenomenological behavior of 
a hydrogen cascade arc. They are: current, pressure and 

tube radius. The heavy particle temperature at the wall 
has been fixed at 1000*K for all results presented. Other 
variables which have been considered are the reflection 
coefficients for electrons ($) and ions (x) > which are 
allowed to vary from a maximum of 0.99 to a minimum of 0.0 
(or somewhat higher depending on the stability of the numeri- 
cal scheme) . Because the influence of the reflection coef- 
ficients on the flow profiles and the integrated quantities 
(li, heat flux, etc.) are relatively minor, they are not 
considered to be important parameters. 

Results typical of those obtained for the influence of 
the reflection coefficients are shown in Figures 4-6 and 4-7 
for the case of I = 200 amps, P = 1 atm and R = 0.005 meters. 
These results illustrate the electron temperature and number 
density profiles for values of $ and y equal to 0.99 and 0.20. 
Notice that 3 and y influence the profiles only in the region 
very close to the wall. As the reflection coefficients 
decrease, more electron-ion pairs are combined into atoms 
at the wall (Figure 4-7), and it is therefore not surprising 
to note a decrease in electron number density at the wall. 

Due to the fact that the heavy particle temperature profile 
is virtually unaffected by 3 and y» the heavy particle heat 
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Figure 4-6. Illustration of Electron Temperature Profile 
Variation as a Function of the Reflection 
Coefficients . 
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Figure 4-7. Illustration of Electron Number Density Profile 
Variation as a Function of the Reflection 
Coefficients . 
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conduction varies by only \% between the two conditions 
shown. Since the radiation level is essentially determined 
by the high temperatures in the arc core, it too is inde- 
pendent of B and x* contrast, electron heat conduction 

and diffusion effects are directly related to sheath param- 
eters and are therefore relatively sensitive to the reflec- 
tion coefficients. However, these terms have been shown to 
be insignificant in the present study since, for all operat- 
ing conditions, their sum accounts for less than 21 of the 
total energy balance and, on the average, for less than 0.5%. 
Finally, it has also been shown that the electric field 
intensity, and therefore the total power input, is also 
virtually independent of the values of $ and x* 
their comparative unimportance, the reflection coefficients 
have therefore been fixed at 0.99 for all further calcula- 
tions . 

The parametric studies are now discussed in three parts. 
These pertain to consideration of the current, the pressure, 
and the tube radius as independent variables. The tempera- 
ture and electron number density profiles calculated for a 
pressure of 1 atm and a radius of 0.005 meters are shown in 
Figures 4-8 and 4-9 for arc currents of 50, 100, 150 and 200 
amps. In Figure 4-8, the upper curve of each pair originating 
from the same centerline temperature is the electron tempera- 
ture and the lower curve is the heavy particle temperature. 

The electron wall temperature varies from 1260°K for 
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Figure 4-8. Electron and Heavy Particle Temperature Versus 
Current at P = 1 atm and R = .005 m. 
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[ = 50 amps to 1930°K for I = 200 amps. The figures reflect 
typical arc behavior. As the operating current increases, 
there is an increase in temperature due to the increased 
energy input. Also, due to increased temperature, there is 
a corresponding increase in electron number density. One 
interesting aspect of Figure 4-8 is that there is little 
difference between the electron and heavy particle tempera- 
ture profiles for all currents, indicating that for these 
conditions thermal equilibrium is closely approached. In 
fact, for all solutions which were obtained, the thermal 
nonequilibrium effect is minor compared with chemical non- 
equilibrium. As mentioned in Chapter 2, the degree of 
thermal nonequilibrium in a gas is dependent upon, among 
other things, the atomic weight of the gas. That is, lighter 
gases exchange kinetic energy with electrons more efficiently 
than heavier gases. It is therefore not surprising to note 
that, for a heavier gas such as argon, Clark [9] reports a 
high degree of thermal nonequilibrium with electron wall 
temperatures "freezing" at values between 6000 and 8000°K. 

Figure 4-10 illustrates the nature of the chemical 
nonequilibrium condition for currents of 50 and 200 amps. 

The solid lines represent the computed nonequilibrium elec- 
tron number density, and the solid lines with circles are 
the electron density calculated from the Saha equation at 
the electron temperature. Severe chemical nonequilibrium 
in the wall region is evident in all results of this study. 
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Figure 4-10. Illustration of Current Dependence of Chemical 
Nonequilibrium. 
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The nonequilibrium is induced by the ambipolar diffusion of 
clectron-ion pairs. In the wall region the net flux of 
electrons to a differential element increases above the rate 
at which electron recombination occurs. The net effect. is 
the elevation of the electron concentrations in the arc 
periphe ry above the values which would exist under equilibrium 
conditions. Note in Figure 4-10 that the radial location 
which marks the onset of chemical nonequilibrium increases 
with increasing current. This is because, at a given radial 
location, increasing the current increases both the electron 
temperature and concentration, which in turn increases the 
rate of the equilibrium restoration processes (in this case 
collisional and radiative recombination). 

In the computer program, the electric field is calculated 
through numerical integration of Equation (3-4). Figure 4-11 
displays the influence of arc current on the electric field 
intensity of the hydrogen cascade arc. It is interesting 
to note that both equilibrium and nonequilibrium models 
predict a local minimum of the electric field. For the non- 
equilibrium curve, the minimum occurs near 130 amps. The 
existence of this minimum may be explained by considering 
the electric field to be equal to the current divided by the 
over-all conductivity of the arc column, as in a solid con- 
ducting rod. An increase in current results in an increased 
degree of ionization and, therefore, an increase in the con- 
ductivity of the arc column. At lower currents, a small 
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Figure 4-11. Equilibrium and Nonequilibrium Calculations of 
Electric Field-Current Characteristic in 
Hydrogen and Argon. 
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increase in current results in a large enough increase in 
the degree of ionization, and hence the conductivity, to 
provide for a decrease in the electric field intensity. At 
higher currents, the degree of ionization is less responsive 
to a current increase and any increase in current increases 
the electric field intensity. As mentioned in Section 4.2, 
the reason for the discrepancy in the equilibrium and non- 
equilibrium pVedictions of electric field intensity is due 
to the inability of the equilibrium model to accurately 
calculate the electrical conductivity outside the core 
region. 

The reason for the hydrogen electric field intensity 
exceeding that for argon is due to the fact that because of 
its greater atomic mass (hence, less efficient transfer of 
elastic energy from electrons to heavy particles) , the elec- 
tron concentration and temperature of argon far exceed that 
for hydrogen outside the core. Therefore, the integral of 
the electrical conductivity is greater and the electric 
field lower in argon. 

Before discussing the heat transfer properties of the 
hydrogen arc, it is important to understand how the heat 
fluxes are computed. In Section 2.2.3 reference was made 
to the method of calculating the electron wall heat conduc- 
tion by consideration of the electrostatic sheath. Without 
derivation, the equations used in calculating the electron 
wall heat conduction are 
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(4-1) 
(4-2) 
(4-3) 

I (4-4) 

where and £ are the ion and electron wall fluxes, 

w ®w 

Cw and c are the ion and electron thermal speeds, J is 

n 0 W 

W W 

the wall current (zero in this study) , A is the electro- 

w 

static wall potential and e is electronic charge. In addi- 
tion, the wall ambipolar diffusion velocity is: 

« % 

Equation (4-4) is used to obtain electron wall heat conduc- 

tion. For atomic hydrogen, the summation over diffusion 
velocities in the heavy particle energy equation (2-9) 
reduces to zero. The heavy particle energy equation, when 
multiplied by the radial coordinate and integrated from 
centerline to wall, then yields, for the heavy particle heat 
conduction , 

f<lr H 

n 'o 




92 


The remaining heat flux quantities are obtained from the 
electron energy equation (2-8). Expressing this equation 
in terms of electron temperature and number density, multi- 
plying by the radial coordinate, and integrating, the equa- 
tion becomes 


kR(n TV.) 

'■ e e amb-^w 


= - R(q ) + ^ + rS dr + rSf dr 

'•^r 27 t J ec J ic 


(4-7) 


o 

The term S. can be rewritten as the sum of a radiation and 

1C 

a three body collision term. When this is done, the three 
body collision term simply becomes the product of the ioniza- 
tion potential and the electron creation rate, Ip^g • Making 
these substitutions and dividing Equation (4-7) by R gives 

1 k(n T V . ) = - (q ),, + ^ - (q^ ’ i f r Rad dr 

2 '■ e e amb-^w '•^r^'^w 27rR Tu w R I 

e n Jq 


1 


' R 



rl n dr 
P e 


(4-8) 


Finally, recognizing that the last term of Equation (4-8) is 
simply the product of Ip and the integral of the electron 
continuity equation (2-7) , the energy equation may be 
rewritten as 


e n J 


R 


‘rh^w 


r Rad dr 


* ^p^^e^amb^w 


(4-9) 
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The Left-hand side of Equation (4-9) is then the total power 
input to the gas and the right-hand side contains the various 
energy loss terms. The first and second terms of the right- 
hand side are the electron and heavy particle heat conductions 
calculated by Equations (4-4) and (4-6). The third term is 
the flux of electron kinetic energy to the wall due to ambi- 
polar diffusion, and the fourth term is the radiation heat 
flux to the wall. Note that, since the cascade arc is 
assumed to be an infinitely long cylinder and the gas is 
assumed to be optically thin, the energy radiated per unit 
length of tube by the gas is equal to the radiation energy 
absorbed per unit length by the tube wall. The fifth term 
is the product of the ionization potential and the ambipolar 
electron flux at the wall and is therefore equivalent to the 
diffusion of ionization energy to the wall. The third and 
fifth terms are calculated with the aid of Equation (4-5) and 
the fourth term by numerical integration. In all of the 
calculations of this study, the combined electron heat con- 
duction and diffusion terms accounted for less than 2% of 
the total energy loss. Due to their minor influence on total 
wall heat transfer, the values of these terms are not plotted 
in succeeding graphs. It cannot be said, however, that 
ambipolar diffusion, electron heat conduction and the diffu- 
sion of ionization energy within the arc do not have a very 
important influence on the temperatures and densities and 
hence, indirectly upon the radiation, heavy particle conduc- 
tion and total power input, of the hydrogen arc. 
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The very nature of the asymptotic region is such that 
the total power input is equal to the heat transfer to the 
wall. Therefore the total wall heat flux is given by the 
term on the left-hand side of Equation (4-9) and is shown 
plotted in Figure 4-12. The reason, then, for the hydrogen 
heat flux exceeding that for argon is the higlier electric 
field of argon. The higher electric field for argon was 
explained earlier in this section as being due to the higher 
electron concentrations and temperatures outside the core 
region. The wall heat flux is primarily composed of radia- 
tion and heavy particle conduction. 

The percentage radiative heat flux is shown in Figure 
4-13. The percentage of heavy particle conduction is approx 
imately 100% minus the percentage radiation. Figure 4-13 
demonstrates that at lower currents, heat transfer from the 
arc is almost completely due to heavy particle conduction. 
However, at higher currents radiation assumes part of the 
burden of transferring heat from the gas. Due to the much 
higher electron concentration and temperature outside the 
core in argon, the radiative recombination rate (which is 
proportional to n^) is higher and a larger portion of the 
argon arc will be radiating significantly. Hence, the radia 
tion percentage is higher for argon than for hydrogen. 

Figures 4-14 and 4-15 reveal the effect of arc pressure 
on the electron temperature and concentration profiles. In 
Figure 4-14, the electron and heavy particle temperatures 
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Figure 4-12. Equilibrium and Nonequilibrium Calculations of 
Total Wall Heat Flux as a Function of Current 
in Hydrogen and Argon, 
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Fieure 4-13. Equilibrium and Nonequilibrium Calculations of 
Radiative Wall Heat Flux Fraction as a Function 
of Current in Hydrogen and Argon. 
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Figure 4-14. Electron Temperature Profiles Versus Pressure 
at I = 200 amps and R = 0.005 m. 



Figure 4-15. Electron Number Density Versus Pressure at 
I = 200 amps and R = O.OOS m. 
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are close enough to be indistinguishable, hence only is 
plotted. The electron wall temperature varies from 1930*K 
for P = 1 atm to 1032°K for P = 10 atm. The decrease in 
centerline temperature with increasing pressure is best 
explained by referring to Figure 4-19. Through this range 
of pressure, the radiation flux increases from ISI at 1 atm 
to about 551 at 10 atm. This increase of radiation transfer 
from the arc core causes the reduction in temperature shown 
in Figure 4-14. Also, as pressure rises, electron-heavy 
particle collision rates increase resulting in a stronger 
tendency toward equilibrium and thus the decrease of wall 
electron temperature and wall number density. The increase 
in electron number density in the core with pressure is a 
result of the increasing atom number density and subsequent 
ionization rate. Figure 4-16 illustrates the nature of the 
chemical nonequilibriura dependence on pressure. The increase 
of electron-heavy particle collision rates with pressure 
suppresses chemical nonequilibrium for a greater distance 
from the centerline. 

Both Clark’s results and those of the present study 
show an increase in electric field intensity with pressure 
as shown in Figure 4-17. Since the current is fixed, this 
implies a decrease in the integrated electrical conductivity 
with pressure. Hence, the positive influence of increasing 
electron number density on the electrical conductivity is 
not as strong as the negative influence of increased elastic 
collision rates (see Equation (C-2)). 
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Figure 4-16. Illustration of Pressure Dependence of Chemical 
Nonequilibrium. 
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Figure 4-17. Nonequilibrium Electric Field Intensity as a 
Function of Pressure for Hydrogen and Argon, 
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Figure 4-18 illustrates the dependence of total wall 
heat flux on pressure both for hydrogen and argon. Since 
the Ohmic dissipation and the heat output must balance in 
the asymptotic region, an increase in the field intensity, 
and hence in power input, results in increased wall heat 
flux for fixed arc current. A substantial portion of the 
increase in wall heat flux with pressure can be attributed 
to the increase in radiative emission, as indicated by 
Figure 4-19. The increased number densities associated 
with higher pressure increase the radiative recombination 
rate and hence the radiative emission. 

Figures 4-20 and 4-21 depict the effect of tube radius 
on the temperature and electron number density profiles. 

The decrease in temperature and electron density with 
increasing radius is due to an accompanying decrease in the 
energy input per unit area as the radius is increased. For 
the 0.005 m radius profile the electron wall temperature is 
lOSO^K, and for the radius of 0.05 m it decreases to 1014'‘K. 
Also as the tube radius is increased, the electron and heavy 
particle temperatures become essentially equal indicating 
that thermal equilibrium is established. This is due to the 
sharp reduction in the strength of those processes which 
promote thermal nonequilibrium, particularly the selective 
transmission of energy to the electrons by the electric 
field, with increasing radius. 
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Figure 4-18. Nonequilibrium Total Wall Heat Flux as a 

Function of Pressure for Hydrogen and Argon. 
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Figure 4-19. Nonequilibrium Radiative Wall Heat Flux Fraction 
as a Function of Pressure for Hydrogen and Argon. 


Te a T^, (*^K)X10 


105 



Figure 4-20. Electron and Heavy Particle Temperature Versus 
Tube Radius at I = 200 amps and P = 1 atm. 
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Figure 4-21. Electron Number Density Versus Tube Radius at 
I = 200 amps and P = 1 atm. 
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The nature o£ the chemical nonequilibrium condition is 
shown in Figure 4-22 for radii of 0.005 and 0.05 meters. 
Ucparture from the chemical equilibrium condition occurs 
nearer to the centerline as the radius increases. This is 
due to the sharp reduction in the local ionization and 
recombination rates (the principal mechanisms for maintain- 
ing chemical equilibrium) that occurs with increasing radius. 

The variation of the electric field intensity with tube 
radius is indicated in Figure 4-23. There are two major 
opposing processes which influence the field intensity as 
the tube radius is increased. With the current held fixed, 
the electric field is inversely proportional to the integrated 
conductivity of the arc cross section. As the radius in- 
creases, the higher cross-sectional area of the gaseous con- 
ductor serves to increase this integrated conductance. How- 
ever, the attendant decrease in the temperature and electron 
density which causes the electrical conductivity to decrease, 
serves to lower the conductance. From Figure 4-23, it is 
apparent that the increasing cross-sectional area effect is 
dominant, causing the electric field to decrease with radius. 

The total wall heat transfer per unit length of tube is 
plotted in Figure 4-24 and reflects the trend of the electric 
field intensity shown in Figure 4-23, since the total heat 
transfer is the product of the current and the electric field. 
These same trends are predicted by Clark [9] for argon. 
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Figure 4-22. Illustration of Tube Radius Dependence of 
Chemical Nonequilibrium. 
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Figure 4-23, 


Nonequilibriura Electric Field Intensity as a 
Function o£ Tube Radius for Hydrogen and Argon. 
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Figure 4-24. Nonequilibrium Total Wall Heat Rate per Unit 
Length as a Function of Tube Radius for 
Hydrogen and Argon. 
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Figure 4-25. Nonequilibrium Radiative Wall Heat Flux 

Fraction as a Function of Tube Radius for 
Hydrogen . 




112 


The percentage o£ radiation heat flux is plotted in 
Figure 4-25, and the results reflect several opposing 
phenomena. As radius increases, the electron and heavy 
particle temperatures decrease, thus diminishing the radia- 
tion emission per unit volume from the arc core. In contrast, 
the increase in the volume of the core region and the decrease 
in the total heat transfer serve to increase the percentage 
of radiation flux with radius. From Figure 4-2S it is 
apparent that the first phenomenon is dominant at lower radii, 
but one or both of the other two dominate at larger radii. 
However, neither of the magnitudes of the radiative heat flux 
or radiative heat transfer per unit tube length exhibit a 
minimum but decrease monotonically as the radius increases. 

4.4 Uncertainty Calculations 
It has been noted at various points that it is necessary 
to use approximate theories and methods to obtain certain 
properties required by the mathematical model. The use of 
approximate theories is necessary due to the lack of accurate 
data regarding transport coefficients for nonequilibrium 
gases and because of the complexity of a rigorous treatment 
of the radiation processes in a real gas. Therefore, it has 
been necessary to utilize a mean free path theory (Appendices 
A, B, C and D) to obtain the transport coefficients and an 
approximate model for treating radiative emission (Section 
4.2). It is therefore important to perturb these properties 
in order to estimate the over-all uncertainty of the predictions. 
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This was done by rerunning computer solutions for which 
pertinent coefficients and properties were altered and com- 
paring the results with those previously discussed. 

The three-body recombination coefficient and the inverse 
ionization coefficient were divided by a factor of two, and 
a solution was obtained for I = 200 amps, P = 1 atm and 
R = .005 m. The factor of two is representative of the 
discrepancies between existing recombination data. Although 
atomic properties can be calculated rather accurately for 
hydrogen, they are usually determined by considering only 
one process at a time (e.g. radiative recombination exclusive 
of collisional modes of recombination). Therefore, even 
though the individual properties may be accurate, they are 
not necessarily additive. 

With the ionization and three-body recombination coeffi- 
cients reduced by a factor of two, the temperatures and 
electron densities in the core were found to be virtually 
unaffected, due to the presence of chemical equilibrium in 
this region. However, the electron temperature and density 
increased somewhat in the wall region. Also, the electric 
field and heat fluxes were essentially unchanged by this 
rather large perturbation. The conclusion then is that 
uncertainties in the ionization and three-body recombination 
coefficients do not introduce significant errors in the 
mathematical model. 
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As discussed in Section 4.2, the radiation model used 
is rather approximate, and to obtain agreement with accurate 
calculations in the equilibrium limit, it was necessary to 
multiply results by a factor of 0.4. For the uncertainty 
calculations it was therefore decided to use the uncorrected 
value of the recombination coefficient. This procedure 
resulted in a decrease of 800°K in the two temperatures and 
a 251 decrease in electron density at the centerline. The 
electric field and total wall heat flux increased by about 
81 and the radiative heat flux increased by 65%. It is 
therefore obvious that an uncertainty in the radiative 
recombination coefficient can have a significant effect on 
the temperature and concentration profiles as well as the 
radiative flux. However, the effect on the total wall heat 
flux and the field intensity remains small. 

The rate of energy transfer from the electron subgas to 
the heavy particles due to elastic collisions is 

another property that needs to be investigated in an uncer- 
tainty analysis. The value of was divided by 2 and a 
solution obtained for I = 100 amps, P * 1 atm and R = .005 ra. 
This resulted in a negligible change in temperatures and 
electron densities in the core of the arc but increased the 
electron temperature in the wall region by 200 to 300®K. 
Electron densities near the wall increased on the order of 
10%. Changes in the electric field and the total and radia- 
tive wall heat flux were negligible. Considering the fact 
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that an unusually large perturbation factor was used to 
analyze the uncertainty associated with S , it appears that 

C w 

more realistic uncertainties in this property (much less than 
a factor of 2) have negligible effect on the accuracy of this 
model. 

In the mean free path derivation of the transport coef- 
ficients, the calculation of the ambipolar diffusion coeffi- 
cients is thought to be the most accurate (see Figure A- 2) . 
Therefore, uncertainties in the diffusion coefficients are 
expected to be small, and in the uncertainty analysis the 
coefficients were divided by a factor of 1.25 for a current 
of 150 amps and tube radius of .005 m. As a result of this, 
temperatures in the core increased by about 3% while remain- 
ing virtually unchanged in the periphery. Electron concen- 
tration increased by 14% at the centerline and decreased by 
about 10% in the wall region. The electric field and total 
wall heat flux were essentially unchanged while the radiative 
flux increased by about 10%. It is therefore felt that 
uncertainties in the ambipolar diffusion coefficient have a 
negligible effect on computed results. 

The uncertainty associated with errors in the thermal 
conductivity was determined by multiplying both the electron 
and heavy particle thermal conductivities by a factor of 
1.25, and a solution was obtained at I = 100 amps, P = 1 atm 
and R = .005 m. The result was a slight decrease in tempera- 
ture and electron density in the core, less than a 5% 
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decrease in temperature near the wall, and a rather large 
(401) decrease in electron density very near the wall. In 
terms of bulk parameters, the electric field and total wall 
heat flux increased by 111 and the radiative flux decreased 
by 201. Uncertainties in the thermal conductivities are 
therefore felt to have a negligible bearing on the results 
of this study. 

Similarly, the electrical conductivity was perturbed by 
a factor of 1.25, and a solution was obtained for I = 200 
amps, P = 1 atm and R - .005 m. This resulted in moderate 
decreases in temperature (less than 51) over the arc cross 
section. However, number densities decreased by about 141 
at the centerline and 401 near the wall. The electric field 
and total heat transfer decreased by about 10% and the radia- 
tion flux by 25% . 

From the previous considerations, it is evident that 
representative uncertainties in the transport coefficients 
and the atomic parameters have, at best, a minor influence 
on the thermochemical profiles and the over-all arc param- 
eters. The quantity which is most sensitive to the uncer- 
tainties is the radiation heat flux, which may be in error 
by as much as 50%. 

4 . 5 Speculations for the Molecular Hydrogen Arc 

The original objective of this study was to mathematically 
model the molecular, rather than the atomic, hydrogen arc. 

To do so involves combining the model of this study with a 
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molecular continuity equation and representations of the 
thermodynamic and transport properties which account for the 
presence of the molecular species. In fact the molecular 
model was derived, and a computer program was developed to 
solve the governing equations of the molecular hydrogen arc 
in the entrance region. A great deal of time was spent in 
devising the computer program, but continued problems with 
numerical instabilities precluded obtaining realistic solu- 
tions. Therefore, exclusive emphasis was placed on study 
of the atomic hydrogen arc. 

It is possible to extrapolate the atomic hydrogen arc 
behavior to that of the molecular arc by considering the 
differences in the properties of the two gases. Note that, 
if the entire arc is in thermochemical equilibrium, the 
molecular number density in the core will be several orders 
of magnitude smaller than that of the electrons or atoms. 
Furthermore, close to the wall, the atom number density will 
be negligible but the atom diffusion velocities will be large. 
From the results of this study, thermochemical equilibrium is 
known to exist in the core of the arc for all operating con- 
ditions and it is reasonable to assume that this condition 
would prevail for the molecular arc. For atomic hydrogen, 
it was also found that the large electron number density 
gradients near the wall promote a departure from chemical 
equilibrium. It is expected that a similar condition would 
be predicted for the molecular hydrogen arc. In fact. 
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relaxation time calculations reveal that characteristic 
electron and atom diffusion times are several orders of 
magnitude smaller than electron and atom recombination 
times. Therefore the electron and atom number densities 
in the wall region of the molecular arc should be many 
orders of magnitude larger than values characteristic of 
chemical equilibrium. Therefore actual molecular hydrogen 
number densities near the wall are expected to be somewhat 
lower than those values characteristic of equilibrium. 

According to the calculations of Yos [21] , Devoto [20] 
and Vanderslice [28], the thermal conductivity of a H-Hj 
mixture is somewhat greater than that of atomic hydrogen for 
the temperatures of interest in this study. As noted in 
Section 4.4, an increase in thermal conductivity tends to 
’’flatten” the temperature profile. Therefore, in the 
molecular arc, the higher heavy particle thermal conductivity 
will tend to flatten the heavy particle temperature profile 
near the wall. 

A factor which will tend to depress the heavy particle 
temperature in the wall region of the molecular arc relates 
to the additional modes of energy storage associated with 
the rotational and vibrational states of the molecule. This 
implies that the molecular gas can store more energy at a 
lower temperature than atomic hydrogen. 

The presence of molecular hydrogen in the wail region 
decreases the efficiency of elastic energy exchange between 
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electrons and heavy particles due to the higher mean heavy 
particle mass. This will act to lower the heavy particle 
temperature and raise the electron temperature. 

As mentioned previously, it is expected that in the 
nonequilibrium molecular arc, the local atom number density 
will be somewhat higher than at equilibrium (though still 
lower than in the atomic hydrogen arc) near the wall. This 
implies an imbalance between the recombination-dissociation 
rates in which recombination dominates. When atomic recom- 
bination occurs to form a hydrogen molecule, the dissociation 
energy is transferred to another heavy particle or an elec- 
tron. Therefore, atomic recombination in the wall region 
tends to elevate the thermal state of the gas. 

Mean free path calculations for the molecular hydrogen 
arc indicate that the ambipolar diffusion coefficients will 
be somewhat lower in the wall region than for the atomic 
hydrogen arc. Since this tends to inhibit the flux of 
electrons and energy from the core to the wall region, the 
effect is to slightly increase electron density and tempera- 
ture in the core and to diminish these quantities near the 
wall . 

Considering all of the above factors, it is felt that 
those which dominate are the higher mean heavy particle 
mass and the greater energy storage capacity of the hydrogen 
molecule. Since both of these effects are influential only 
near the wall, it is felt that the nonequilibrium molecular 



120 


hydrogen arc will exhibit lower heavy particle temperatures 
and higher electron number densities and temperatures in the 
wall region. The more minor effect of lower ambipolar dif- 
fusion velocities near the wall may also promote a small 
increase in both temperature and electron density at the 
centerline . 

Assuming the above to be true, the increased electrical 
conductivity in the wall region, brought about by the higher 
electron temperature and number density, will result in a 
larger integrated conductivity and, hence, lower electric 
field in the nonequilibrium molecular hydrogen arc as com- 
pared to the atomic hydrogen arc. Furthermore, since the 
total energy input to the arc, which is the product of 
electric field and current, must be equal to the total wall 
heat flux in the asymptotic region, the molecular arc should 
exhibit somewhat lower wall heat fluxes due to the smaller 
electric field. Finally, the radiative emission of the 
molecular arc may be slightly augmented due to small in- 
creases in the centerline electron number density and tem- 
perature and also due to molecular radiation effects in the 
infrared spectrum. 

In summary, the nonequilibrium molecular hydrogen arc 
is expected to exhibit higher electron number density and 
temperature and lower heavy particle temperature near the 
wall, slightly higher electron number density and temperature 
in the core, somewhat smaller electric field and total wall 
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heat flux, and slightly larger radiative emission as compared 
to the nonequilibrium atomic hydrogen arc. However, as far^ 
as the bulk quantities are concerned, such as electric field 
and heat fluxes, it is expected that the solutions to the 
nonequilibrium atomic hydrogen arc provide accuracy suffici- 
ent for engineering purposes. 
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CHAPTER 5 SUMMARY AND CONCLUSIONS 

In this chapter, the key results o£ the study are 
summarized. For ease of reference, the major results and 
conclusions are enumerated more or less in the same order 
that they appear in Chapter 4. 

1. Fairly good agreement is obtained between the theoreti 
cal results of this study and the experimental tem- 
perature profiles of Wiese [27] and Steinberger [12] 
for cascade arcs of .0015 and .0010 m radius, respec- 
tively, and a current of 40 amps at atmospheric 
pressure . 

2. Comparison of nonequilibrium predictions with the 
experimental electric field intensity and total heat 
flux data of Wiese [27], Morris [10] and Maecker [11] 
for their small radius arcs gives poor agreement at 
low currents with some improvement at higher currents. 
The reason for this disparity relates to convergence 
errors in the second order finite difference operator 
which are particularly severe at small tube radii. 

3. Differences in the nonequilibrium and equilibrium [6] 
predictions for the hydrogen arc are readily attribut- 
able to an underprediction of the electrical 
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conductivity by the equilibrium model in the wall 
region. This results in the equilibrium model pre- 
dicting high centerline temperatures and low tempera- 
tures in the wall region. 

4. Use of the equilibrium theory results in a slight 
overprediction of the field intensity and the total 
wall heat flux. 

5. Comparison of nonequilibrium and equilibrium radiative 
heat fluxes gives good agreement at higher currents, 
but the equilibrium results are considerably higher at 
lower currents. 

6. The value of the electron and ion reflection coeffici- 
ents has very little influence on the temperatures 
and electron density, except immediately adjacent to 
the wall. Electric field intensity and the total and 
radiative wall heat fluxes were unaffected by variation 
in these parameters. The reflection coefficients do 
have a strong influence on the electron heat conduction 
and the diffusion of kinetic and ionization energy at 
the wall, but the sum of these three terms accounts for 
less than 21 of the wall heat transfer for all parametric 
s tudies . 

7. When the arc current is increased with all other 
operating parameters fixed, both temperature profiles 
and electron concentration profile increase over the 
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entire arc cross section. The radial location of the 
onset of thermal and chemical nonequilibrium increases 
with increasing current. The electric field exhibits 
an interesting phenomenon by first decreasing, passing 
through a minimum, and then increasing with current. 

The total, radiative, and heavy particle conduction 
wall heat fluxes increase with increasing current. 

8. When the pressure is increased and all other operating 
parameters fixed, the centerline and wall temperatures 
decrease significantly, the electron concentration 
increases at the centerline and decreases at the wall 
and the radial location of the onset of thermal and 
chemical nonequilibrium increases. Both the electric 
field intensity and the total wall heat flux increase 
with pressure, and the radiative heat flux percentage 
increases from 151 to 55% as pressure increases from 
one to ten atmospheres. Heavy particle conduction 
actually decreases slightly as the pressure is increased. 

9. When the tube radius is increased with all other 
operating parameters fixed, both temperature profiles 
and the electron concentration profile decrease over 
the entire arc cross section and the normalized radial 
location marking the onset of thermal and chemical 
nonequilibrium decreases. For tube radii I greater than 
.02 m and I = 200 amps, P = 1 atm, thermal equilibrium 
is established over the entire arc cross section. The 
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electric field and the total and heavy particle con- 
duction wall heat transfer all decrease rapidly with 
increasing tube radius. Although the radiative heat 
transfer percentage passes through a minimum, the 
radiation heat transfer decreases rapidly as tube 
radius increases. 

10. Although chemical nonequilibriura is never significant 
near the tube centerline, it always exists in the wall 
region. 

11. Thermal nonequilibrium is relatively unimportant for 
the operating conditions considered in this study. 

The highest electron wall temperature is 1930°K which 
occurs for I = 200 amps, P = 1 atm and R = .005 m. 

12. In general, comparison with the nonequilibrium calcula- 

tions for argon [9] reveals the same basic trends, 
with the following exceptions: thermal nonequilibriura 

is a much more pronounced effect in argon and, in argon, 
the diffusion of thermal and ionization energy and 
electron heat conduction are not negligible contribu- 
tions to the total heat flux. 

13. An uncertainty analysis revealed that reasonable errors 
associated with the utilization of ionization and 
three-body recombination coefficients in the inelastic 
collision model had an insignificant effect on the 
solutions. However, it was discovered that uncertainties 
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in the radiative recombination coefficient, although 
having negligible effect on electric field intensity 
and total heat flux, could result in errors in the 
radiative heat flux as high as 651. Also, uncertain- 
ties in ambipolar diffusion coefficients had an insig- 
nificant effect on the predictions. Reasonable 
inaccuracies in the atomic and transport properties 
used in this study may result in uncertainties of the 
order of 1 % in centerline and wall temperatures, 25S 
and 40% in centerline and wall electron concentration, 
11% in the electric field intensity and total wall heat 
flux and 50% in the radiation heat flux. 

14. Qualitative considerations for the nonequilibrium 

molecular hydrogen arc suggest that atom and electron 
concentrations should be much higher and molecular 
concentrations somewhat lower than those characteristic 
of equilibrium in the wall region. Furthermore, in 
comparison with the nonequilibrium atomic hydrogen arc, 
it is felt that the molecular hydrogen arc will exhibit 
lower heavy particle temperatures and higher electron 
number densities and temperatures in the wall region. 
Also, electron concentration and temperature are 
expected to be slightly higher in the core. Finally, 
it is expected that, although electric field intensity 
and total wall heat flux will decrease and radiative 
heat flux increase somewhat in the nonequilibrium 
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molecular arc, the predictions available for these 
quantities from the atomic hydrogen solutions provide 
accuracy sufficient for engineering purposes. 
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APPENDIX A 

AMBI POLAR DIFFUSION IN HYDROGEN BY MEAN FREE PATH THEORY 

The purpose o£ this appendix is to describe the theory 
behind calculation o£ transport coe££icients by mean £ree 
path methods and to derive an expression £or the ambipolar 
di£fusion velocity (see assumption 13 of Section 2.1). 

Consider an imaginary plane fixed with respect to the 
walls of a container in which gaseous, atomic hydrogen 
exists (see Figure A-1). It is desired to calculate the 
thermal flux of particles passing through this plane in the 
absence of gradients of temperature and concentration. The 
number of particles with velocity vectors lying in the direc- 
tion of a cone of solid angle dw per unit volume of gas is 
n ^ where n is the number density of particles. If 0 is 
the angle between the normal to the plane and the line 
through the center of dm, the component of velocity of the 
particles (with velocity vectors in dm) normal to the plane 
is c COS0 where c is the mean thermal speed of the particles. 
Therefore, the flux of particles through the plane with 
velocity vectors lying in dco is n ^ c cos0. Integrating 
over the hemisphere gives 

ir/2 1 _ 

sine cose ded<() = ^ nc (A-1) 

■ 0 '0 

which is the thermal flux of particles through the plane. 


nc 

TtT 
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in the presence o£ property gradients, the flux of 
particles through the plane is actually the value of ^ nc 
evaluated one mean free path below (or above for negative 
flux) the plane. This is because, on the average, each 
particle arriving at the plane from below does not suffer 
any collisions in its journey to the plane from one mean 
free path away. Using a Taylor series expansion to the 
first order, the flux of particles through the plane from 
below (positive flux) and above (negative flux) is then, 
respectively , 


positive flux 
negative flux 
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(A- 2) 
(A- 3) 


where X is the mean free path length. Therefore, the net 
flux of particles in the positive direction is 
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The free diffusion velocity of particles in the positive 
direction is then given simply as 




3n 

3r 



In a system of three species (electrons, ions, atoms), 
the above expression may be used to obtain the species diffu- 
sion velocities relative to a reference frame fixed in the 


container. 
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Utilizing the assumption that electron and ion number 
densities are equal (Section 2.1), the three diffusion equa- 
tions are 


rD' = 
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However, the diffusion velocity of each species is 
influenced by the presence of the other species, an effect 
which is not accounted for in the above three equations. 

In order to remedy this situation, a fourth equation is 
introduced to force the mass average of the diffusion 
velocities to be zero: 


Pe 
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Ph 
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The constant a is introduced in such a manner that each of 
the terms in parentheses gives the correct diffusion velocity 
for the respective species. After simplification and use of 
number densities, the above expression yields 


a = - 


n^+nni im^ 


"*e , .D ’ 

=- "e'^e 


^ I/D' 

e H+ 


* Vh' 


The correct diffusion velocities are now written with- 
out the prime superscript and satisfy the following relation- 
ship : 
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(A-5) 


Therefore, the diffusion velocities are calculated as 
follows: 


a • ' 

^u+ 


One of the assumptions of Section 2.1 states that the 
electrons and ions will diffuse, effectively, toward the 
wall from the centerline of the arc in pairs. In other 
words the electron and ion diffusion velocities are equal. 
Since the present derivation has not incorporated the 
assumption of ambipolar diffusion, the above diffusion 
velocities for electrons and ions are not identically equal. 
To incorporate the phenomena of macroscopic electric 
neutrality and ambipolar diffusion, the results of Holt and 
Haskell [13] are used. They give the following equations 
for the flux of electrons and ions in the presence of an 
electric field (from a moment of the Boltzmann equation) : 

Te = = - D^vn^ - n^u^E (A- 6) 

''h* ‘ CA-7) 

where the D's are diffusion coefficients, E is the externally 
applied electric field, and p denotes the charged particle 
mobility defined as 

l®l j |e| 

y = li ^ = 
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In the above equations, is the sum o£ collision frequen 

cies involving electrons and all heavy particles and 
is the sum of collision frequencies involving ions and all 
heavy particle species. By eliminating the electric field 
through simultaneous solution of Equations (A-6) and (A-7), 
the following expression is obtained: 
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'H+ 
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Holt and Haskell define the term in parentheses as the 
ambipolar diffusion coefficient. Hence, the left hand side 
of the equation is the ambipolar diffusion velocity: 
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(A-8) 


Note that the above definition of the ambipolar diffu- 
sion velocity accounts only for ordinary diffusion involving 
the electron (ion) number density gradient. Since the true 
electron and ion diffusion velocities previously defined 
involve other gradients as well, it is assumed that the 
above relationship between the ambipolar and species diffu- 
sion velocities is valid for other modes of ambipolar diffu- 
sion (i.e. thermal diffusion and pressure diffusion). 

Now, the electron and ion diffusion velocities may be 
replaced with the ambipolar diffusion velocity and Equation 
(A-5) is no longer identically satisfied. However, with 
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the ambipolar diffusion velocity having been established, 
the atom diffusion velocity may now be computed to insure 
satisfaction of Equation CA-5). Hence, 

P , + p , + p„vS = 0 (A- 

^e amb amb H 

where is now the correct atom diffusion velocity in a 
gas where ambipolar diffusion of electrons and ions occurs. 
From Equation (A-9) it is obvious that the expression for 
the atom diffusion velocity is 

H nj^ amb 

The expression for the electron-ion ambipolar diffusion 
velocity is given by 

n 3n„ 3T^ 3T, 
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The mean free paths (A) in the above expressions are dis- 
cussed in Appendix D and the mean thermal speed (c^) o£ 
species i is calculated from 


c . 

1 



1/ 2 


In order to compare the previous results with those 
obtained from equilibrium kinetic theory, the coefficients 
must be evaluated in the limit of thermochemical equilibrium. 
To facilitate this comparison the following expressions are 
introduced : 


where P and P 
e 

and X and x„ 
e n 

respectively . 


n. 


p = — h — P = kn T = x^P 

e 2n +nu e e e 

e n 


CA-11) 


P„ = ■= — P = kn,,T, = x„P 

H 2n +n,, H h H 


(A-12) 


‘e “H 

pj are the electron and atom partial pressures 
are the electron and atom mole fractions, 
Hence , 


n. = 


P 

FT 




Fr“ 




Substituting the above expressions into Equation (A-10) , we 


obtain 


V 


D 

amb 


p 

^ , f P 1 


ar 


3Xh 

3r 


fP^e ) 

A A 1 ^ 
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'^h 
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9r 
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The last two terms (in brackets) are associated with the 
thermal diffusion coefficient and for the sake of comparison 
with results generated from the rigorous kinetic theory [20], 
they will be dropped and only ordinary diffusion considered. 
Substituting once more from Equations (A- 11) and (A-12) and 
algebraically manipulating, we obtain 


D 

^Lb = A 


_ e , , 

1 Xg 3r ^2 9r 


= Ai ^ - 2Ao •— — 

^ X 3r ^ X,, 3r 

e H 


Ai 




- 2A, 




n 


H 


3x^ 

9r 


r r^vHjfy- cAi - 2 a.)J 


9r 


(A-13) 


Note that in the second step 


9xh 


= - 2 


3x, 


9r “ 9r * 

In the equilibrium limit for equal electron and heavy 

particle temperatures and electron number densities calcu- 
lated from the Saha equation, the term in brackets in Equa- 
tion (A-13) becomes the equilibrium ambipolar diffusion 
coefficient (D^jj^j^) [23]. 

Figure A-2 is a comparison between results obtained 
from rigorous kinetic theory [20] and from the mean free 
path calculations of this study evaluated in the equilibrium 
limit. Although agreement between the two theories is good 
over the entire range of temperature, a multiplicative factor 
of 1.2 was introduced to force excellent agreement of the 
mean free path and rigorous kinetic theory results. 
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Figure A- 2. Comparison of Equilibrium Ambipolar Diffusion 
Coefficient from Mean Free Path Calculation 
and Rigorous Kinetic Theory. 
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APPENDIX B 

HEAT FLUX AND SHEAR STRESS IN HYDROGEN 
BY MEAN FREE PATH THEORY 


In this appendix, mean- free-path expressions are 
derived for the electron and heavy particle heat fluxes and 
the total shear stress. 

The flux of particles through an imaginary, infinitesi- 
mal plane (see Figure A-1) in the positive and negative 
directions, is derived in Appendix A, and the results given 
by Equations (A-2) and (A-3). Assume that some property, 

4>, is associated with the particles passing through the 
plane and that its value is a function of the radial 
coordinate. The positive and negative fluxes of that 
property (kinetic energy, momentum, etc.) through the plane 
can then be written from Equations (A-2) and (A-3) as 


positive flux of 4* 



9n , 
3r ^ 




. 

9r 


negative flux of ij* 


1 


n + 



c + 


9c 

9r 


+ 


9 (|) 

9r 


X 


The net flux of the property, 4), in the positive radial 
direction is then obtained by subtracting the negative flux 
from the positive flux 
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net flux = - 


I c<j)X ^ 


- j n<})X 


3c 

3r 


I 


9(J> 

3r 


(B-1) 


where all higher order terms in X have been excluded, n is 
the particle number density and c is the mean kinetic speed 
of the particles. 

To obtain the heat flux contribution of a given species 

in the gas, let 4) be the mean kinetic energy of that species 

For electrons and other non-molecular components, the mean 

5 

kinetic energy per particle is y kT. Making this substitu- 

f 8kTl ^ ^ ^ 

tion for cf) and substituting the expression for c - ^nT”! 

in the derivative of c with respect to r in Equation (B-1) 
gives the following for the radial heat flux of species i 

3T. 




I W 


_ 3n . 


where the first and second terms relate to Fourier conduction 
and thermal diffusion, respectively. In the electron energy 
equation, the electron heat flux is obtained from the above 
equation simply by replacing the i subscript with the sub- 
script e. In the heavy particle energy equation, the heavy 
particle heat flux is obtained by summing the above expres- 
sion over the ions and atoms (i = H'*’ and i = H) . 

In order to obtain the Fourier conduction for hydrogen 
in the equilibrium limit, the electron and heavy particle 
temperatures are equated, the electron number density is 
related to the temperature and atom number density through 
the Saha equation and the Fourier terms of the individual 
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heat fluxes are suninied. Therefore, the equilibrium Fourier 
conductivity of hydrogen from mean free path theory is 
obtained from the expression 


(q ) 

'■^r^ eq 




+ X .n^c„ + 
e e H 




3T 

9r 


where n is given by the Saha equation and the term in 
square brackets is the equilibrium thermal conductivity. 

Figure B-1 provides plots of the above result (labeled 
"mfp theory”]) and that obtained from the rigorous kinetic 
theory calculation for atomic hydrogen thermal conductivity 
[20], Note that if a weighting factor of 3 is multiplied 
by the mean free path expression, very good agreement is 
obtained with the rigorous calculations in the equilibrium 
limit. For this reason, a weighting factor is introduced 
into the nonequilibrium heat flux expressions used in the 
computer programs of this work. 

To obtain the contribution of a given gaseous species 
to the total axial shear stress acting on the radial plane, 
let (j> in Equation (B-I) be the average axial momentum of a 
species i particle, <(> = m^ U, where U is the mass average 
axial velocity. Since shear stress only involves velocity 
gradients, the first two terms of Equation (B-1) are dropped 
giving 


where 

axial 


- -^rzi is 
momentum 


- T 


rz . 
1 


1 n-c. X-m. 

2 1111 3r 


the shear stress (or the viscous flux of 
in the radial direction) and the negative 



(K)„, (WATTS/METER -“K) 
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Figure B-1. Comparison of Equilibrium Thermal Conductivity 
from Mean Free Path Calculation and Rigorous 
Kinetic Theory. 
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sign is introduced on the left side of the equation to 
comply with the usual sign convention. Therefore, the term 
in parentheses is the contribution of the i^^ species to the 
total gas viscosity, which is obtained by summing over all 
species. Therefore 

y = j[mgngCgXg + Vh^H^h] 

where y is the total gas viscosity. In order to obtain the 
equilibrium viscosity, the above expression must be evaluated 
for equal electron and heavy particle temperatures and the 
electron number density must be calculated from the Saha 
equation. ^ 

Figure B-2 is a plot of the equilibrium limit of the 
mean free path calculation of atomic hydrogen viscosity com- 
pared with results of a rigorous kinetic theory calculation 
[20]. Note that in the case of viscosity, a weighting 
factor of about 1.8 would bring the mean free path calcula- 
tions into relatively good agreement with the rigorous 
results. Therefore, a weighting factor of 1.8 is introduced 
into the nonequilibrium shear stress used in the computer 
programs of this work. 

Expressions for the mean free paths used to calculate 
thermal conductivity and viscosity for each gaseous species 
are given at the end of Appendix D. Note that, for these two 
transport phenomena, the value of Jl is 2 for the collision 
frequencies which appear in the mean free path definitions. 



(/Li.)gq, (KILOGRAM /METER- SECOND) X 10* 
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Figure B-2. Comparison of Equilibrium Viscosity from Mean 
Free Path Calculation and Rigorous Kinetic 
Theory . 
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APPENDIX C 

ELECTRICAL CONDUCTIVITY 


To obtain a nonequilibriuiu expression for the electri- 
cal conductivity, it is assumed that the electric field is 
constant and uniform in the axial direction of the cascade 
arc. The force that the electric field exerts on a test 
electron in the axial direction is assumed to be balanced 
by the average rate of electron momentum loss in the axial 
direction (Newton’s Second Law) due to elastic collisions 
with heavy particles. In other words, although each 
individual electron is continually being accelerated by the 
electric field, the effect of elastic collisions with heavy 
particles is to prevent the electron subgas from accelerat- 
ing on a macroscopic scale. A reference frame translating 
with the mass average gas velocity in the axial direction 
is used and the heavy particles are assumed to be stationary 
relative to the energetic electrons. 

Since a heavy particle remains essentially stationary 
while an electron collides with it, the electron momentum 
loss in its original direction of motion is, on the average, 
m c (l-cosy) , where \ is the angle at which the electron is 
deflected away from its original path. 
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The density of electrons with velocities in the vector 
range c , dc is given in Appendix D by Equation (D-1) . The 
flux of such electrons incident upon a stationary heavy 
particle species per unit volume is obtained by multiplying 
the above referenced expression by c^ and nj^ (heavy particle 
number density) to obtain 

The rate at which these electrons pass through impact 
parameter ring segments of area b db db centered about each 
heavy particle (see Appendix D) is 

“ nj^n^f (Cg)Cga(x, Cg)sin(x)dx d<^ dc^ (C-1) 

where a(x, c ) is the differential scattering cross section. 
The momentum loss rate per electron (with velocity in 
Cg, dCg) is obtained by dividing the above equation by n^> 
multiplying by the momentum loss per collision, m^c^ (l-cosy) > 
and integrating over the angles x siid <(). This gives 

nhf(Cg)CedCe I a(x, Cg)sin(x) [mgCg(l-cosx) ]dxd(J) 

where Qgj^CCg) is the electron-heavy particle total cross 
section for diffusion. Integrating over all velocity ranges 
gives the total momentum loss rate of an electron to the 
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heavy particle species 



The momentum loss rate in the axial direction is then 


■"e"h 


z 


The above expression is simplified by replacing Cg , 

z 

the axial component of electron velocity, by the average 
axial component of electron velocity, which is simply the 
electron drift velocity, U^. Therefore, can be removed 
from the integral which results in 

'"eVd|^('"e)<=e ““S 


If now the vector differential dc^ is expanded into the 

scalar differential in spherical coordinates, integration 

is carried out over the angular spherical coordinates , and 

2 

the substitution y = made (see Appendix D) , the 

e 

following expression is obtained: 


”'e"h‘’d 

8 

kTgj 

1/2 ^ 

Zirm 


. 

e 

k J 

^0 J 


= 8 m n, U, (T^) 

e h d eh e-' 


where flg^’**^(Te) is the gas kinetic collision integral for 
electron-heavy particle collisions. As mentioned in Appendix 
D, the collision integral is more commonly avail- 

able and its use will introduce a factor of 0.5 in the above 
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expression. Therefore, the momentum loss rate per electron 
now becomes 

^ "^e ^h ^eh' ^ " "'e '^eh 

where is the electron-heavy particle collision frequency 

for diffusion. 

It was mentioned previously that the electron momentum 
loss rate is balanced by the force exerted by the electric 
field on the electron. This is expressed mathematically as 

eE^ = 

Therefore, U, = - — 

d 


where e is the electronic charge and is the axial elec- 
tric field component. Note that the specific heavy particle 
species (ions and atoms) have been introduced in the expres- 
sion. 

The current density, J, can be defined as 


U 


d 


!e! 

e n 


e 


Therefore , 
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Now that there are two formulas for U^> 
expression can be obtained for the electrical 






an explicit 
conductivity 

(C-2) 


This is the same expression that Holt and Haskell [13] 
obtain by taking moments of the Boltzmann equation. 

By equating the electron and heavy particle tempera- 
tures and evaluating the electron number density from the 
Saha equation. Equation (C-2) can be expressed in the 
equilibrium limit and compared with accurate kinetic theory 
calculations of electrical conductivity. Figure C-1 shows 
the equilibrium limit of Equation (C-2) compared with 
Devoto's [20] calculations of electrical conductivity. 
Agreement between the two curves is good over the entire 
range and, consequently, no correction factors are needed 
to adjust Equation (C-2). The collision frequencies 
appearing in this equation are discussed at length in 
Appendix D. 



(I /OHM -METER) 


152 



Figure C-1. Comparison of Equilibrium Electrical Conductivity 
from Equation (C-2) and Rigorous Kinetic Theory. 
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APPENDIX D 

MEAN FREE PATH DERIVATION 


The mean free path length is an important property in 
the calculation of transport coefficients by the mean free 
path theory. Once the appropriate collision frequencies are 
known it is a simple matter to obtain the mean free path 
length of a given particle, as it is simply the ratio of 
the thermal speed of the particle to the sum of the collision 
frequencies. Therefore, expressions for collision frequency 
will first be derived. 

The Maxwellian velocity distribution function is defined 

such that the density of particles of some species A with 

velocity vectors that lie in the interval (c. , c. + dc. ) , 

..., (c^ , c^ + dc^ ) or c^, dc^, is given by 
z z "z 

Consider now a single particle of another species B, 
having a different velocity vector interval. The vector 
interval of a B particle is Cg, dCg and the magnitude of 
the velocity difference of A and B particles is 
I^B - c^l = Therefore, the flux of A particles 

approaching a single B particle per unit time and area, 
where I | is the magnitude of the velocity of A particles 
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relative to the B particle, is 

The total flux of A particles approaching all B particles 
per unit time then becomes 

®AB I 

Imagine now an annular ring (impact parameter ring) 
segment of radius b, width db and arc d((i centered about each 
B particle. The number of A particles with center of mass 
passing through this impact parameter ring segment of area 
bdbd(|> for all of the B particles per unit volume and time 
is 

nAf(cA)dcA|gAB|nB£(cB)<icB b db d<^ 

However, the quantity bdb is equal, by definition, to 

sinydy where x is the deflection angle of the 
A particle . traj ectory relative to the B particle and 
o ( I gAB I > x) (the differential elastic cross section) is the 
probability of an A particle of relative velocity I g^B I 
being deflected through the angle x* Therefore, the number 
of A particles being deflected into the interval 
(X, X * dx) , ((}>, 4> * d(J)) by B particles per unit time and 
volume is 

^a"b^^^A^^^^B^ ISabI^HIabI » x)sinxdxd(|)dcAdCg 
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The rate at which one A particle "collides" with B particles 
per unit volume resulting in a deflection x>X‘*’dx passing 

through the azimuthal increment 4>,(t>+d(J) is then 

nBfCcA)f(cB) >X) sinxdxd(j)dcAdcB 

If a collision were defined as an interaction between 
an A and B particle which resulted in any nonzero deflection 
X, then the collision rate would be obtained by integrating 
the above expression. Since x is zero, theoretically, only 
for a distance of approach, b, which is infinite, integration 
of this expression might give an infinite result. 

For the purpose of transport property calculations, it 
is desired to include a weighting factor in the above 
expression so that the more important collisions are weighted 
more heavily than the rest. It will turn out that the trans- 
port coefficients are inversely proportional to the collision 
frequency. Therefore, those collisions which tend to promote 
the transport process should be given the smallest weighting 
factor and vice versa. In the case of diffusion, small 
angular deflections promote the process most effectively and 
so a weighting factor of (l-cosy) is appropriate. For 
thermal conduction and momentum transfer (viscosity) , small 
and large deflections (as opposed to 90° deflections) are 
most effective and a weighting factor of (1-cos^x) is appro- 
priate. Other weighting factors could have been selected 
but from rigorous kinetic theory, the above are known to be 
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correct. The normalization factors for these weighting 
factors are 


NF, = 




1 fTT 

(l-cosx)dx 

h 

fn 

I dx 

D 

fTT 

J (l-cos^X)dx 


= 2 


It is the convention to designate the parameter Z as having 
the value of 1 or 2 to denote the diffusion process or the 
thermal conduction and momentum transfer process, respective- 
ly. Therefore, utilizing the normalization factors, the 
weighting factor for diffusion becomes Jl(l-cos^x) > ^ ~ 
and that for conduction and viscosity is il(l-cos^X)> Z = 2. 

The collision rate between one A particle and all B 
particles (with velocity increments c . , dc. and c„ , dc^, 

A A B B 

respectively) for all angles x snd (p becomes 


-»■ ^ ->■ 


<JC|g^Bl,X)a-cos“x)sinxdxd* = QAB(liABl) 


0 0 

where Jl = 1 for diffusion and £ = 2 for thermal conductivity 
and viscosity. By definition, though, 

• 2 7T i-TT 
0 "^0 

which is termed the total cross section. 

Finally, integrating over both A and B particle velocity 
distributions yields the total collision frequency of an A 
particle with the B particles. 
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Afl 


= n 


B 


f(c^)£(c 


B 


) I «AB I ' I «AB I > 


(D-2) 


There are two special cases of the above equation which 
are important. 

^ • — il. ^ electron (e) and B a heavy particle (h) 

Equation (D-2) for this case becomes 


V 


SL 

eh 


n. 


f(<^e)ft‘=h)llehUQeh(8eh’‘*V^h 


Since the electron thermal speed is much larger than that 
of the heavy particles, the above equation can be simplified 
to 

'’eh ' 


According to Vincenti and Kruger [3] the vector differ- 
entials (dcj^ and dc^) can be transformed into scalar differ- 
entials in three dimensional spherical coordinates which 
yield 

2 

f Cj^dcj^sined0d<f» 

0 '0 ' 0 



Hc^) HQl^(ic^)cldc^sinedQd<i> 

0 0 0 

By substituting the Maxwellian distributions and carry- 
ing out integrations wherever possible, the following result 
is obtained. 


V 


eh 


4Trnj^£ 


m 


3/2 


meCg 

■irrr 
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Finally, by making the substitution 


m 

e e 




the above equation assumes the form 

I 


^eh ~ ^*^h^ [zTrra 


r kT^.i/2 




Qeh^S^Y dy 


CD-3) 


II. A and B are particles of nearly equal mass 

In Chapter 2 of their book, Vincenti and Kruger [3] 
prove that if is the velocity of the center of mass of 

the two particles (A and B) , then 


and 



1 2 _ 1 
7 ’"b'^B " I 


Cma + 2 ^®AB 


where 


y 


"^a»"b 
"‘a "‘b 


and 


XM 


"‘a ^ 


mr 




m 


B^B^ 


Assuming the temperature of A and B particles to be equal, 
substituting the above relations into Equation (D-2), and 
carrying out all possible integrations, the following 
expression is obtained for the collision frequency of an A 
particle with a B particle 


yg 


'ab TriTr 


y IV' 


AB 


47T 




0 
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Upon making the substitution 


Y' 


“7W 


the above equation reduces to 


^AB 


kT 


1.2itpJ 


1/2 


I y -V/ 3 n ^ 




CD-4) 


Hirschfelder , et al. [19], give the following definition 
of the kinetic theory collision integral 

Utilizing the above definition for the collision integral 
in the two expressions just obtained for the e-h and A-B 
collision frequencies gives 


^eh 




'-Ib = «"b‘‘«ab"(’') 


where 1=1 indicates diffusion and 1=2 indicates conduc- 
tion and viscosity. 

The significance of the last substitution is that the 
collision (or omega) integrals are important properties for 
the gas kineticist; hence, they are readily available (as 
functions of temperature only) in the literature. However, 
collision integrals are usually only tabulated for the case 
when I = s. This involves the introduction of the factor 
Y in Equations (D-3) and (D-4) . From calculations of 
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(which are readily available in 
analytical form) it is evident that the introduction of the 
factor bas the effect of increasing the collision fre- 

quency by a factor of 2 for s = 1 and by a factor of 4 for 
s = 2. Therefore, the collision frequencies must be divided 


by a compensating factor of 2^, and the final form of the 
collision frequencies for the two special cases is (note 


that = y for Jl 

= $ 

- 1, 2) 




£ 

^eh 

4^h 


II 

III 

* 1 

, 2 

£ 

^AB 

4nB 


£ = s 

= 1 

, 2 

In this study. 

the 

collision 

integrals 

for 

the electron 


atom, ion-atom and atom-atom elastic collisions were obtained 
from the work of Devoto [20] and Yos [21] and the collision 
integrals for charged particle collisions (electron-electron, 
electron- ion , ion-ion) were obtained analytically through 
use of Rutherford's formula for the differential scattering 
cross section in a coulomb field. 

Finally, it remains to define the mean free paths for 
the species of the gas. As mentioned earlier, this is 
simply the ratio of the mean thermal speed of the particle 
(c^) to the sum of the appropriate collision frequencies. 


Therefore 
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A" = 


“T~ — I — rnr 


eH 


5, 


A,,+ 


Ch+ 




^^'1 + 14 + '*’ '^ U * 


r 

H+H 


> = 


^HH-^ ^ 


V, 


HH 


where 


^i = 


SkT^. 1/2 
~irm7^ 


The superscript i on the mean free path definitions 
indicates that the mean free path of a certain species for 
diffusion is not necessarily the same as for thermal conduc- 
tion and viscosity. Also, the mean free paths of ions and 
atoms do not account for collisions of these heavy particles 
with electrons. The reason for this is that when a heavy 
particle collides with an electron, the heavy particle is 
virtually uninfluenced by that collision. This is further 
substantiated by the ’’persistence of velocity" theory 
described by Jeans [22]. 
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APPENDIX E 

ELECTRON 'ATOM BREMSSTRAHLUNG 


The purpose of this appendix is to derive an expression 
for the rate of electron energy loss due to the acceleration 
of electrons in the potential field of a neutral atom. This 
loss of energy due to acceleration of charged electrons will 
herein be referred to as electron- atom Bremsstrahlung. The 
derivation of electron-ion Bremsstrahlung is given by Clark 

[9]. 


Larmor's formula for the rate of radiative emission 
from an accelerating charged particle is [26] 


^^r _ 2fe^a^] 

wr 3[ cM 


(E-1) 


To obtain the energy radiated during a single interaction 
between an electron and an atom, Equation (E-1) can be 
evaluated as an average during the collision and then multi- 
plied by the time of significant interaction. This time of 
interaction is the distance traveled by the electron during 
significant interaction divided by its speed. The interac- 
tion distance is somewhat arbitrarily estimated to be twice 
the atomic (or molecular) diameter and the interaction time 

Z Cl p OTT) 

is therefore expressed as — — , where Cg is the electron 
speed. The energy radiated in a single collision is then 
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dP^ 

Ht“ 


At = 


2 

3 


fe^a^ 


^^atom 


The average acceleration o£ the electron during colli- 
sion is given by the vector difference in electron velocity 
before and after collision divided by the interaction time. 
Since the atom is much heavier and slower it can be con- 
sidered stationary relative to the electron; hence the 
electron speeds before and after collision are approximately 
equal. This implies that the absolute value of the velocity 
vector difference is given from the following diagram as 


z 




2sin 




where \c^' \ = |c^| 

and X is the angle by which the electron is deflected after 
collision with the atom. Finally, the acceleration is given 
by the following formula: 

2cf(l-cosx) 


dP o 2 c" 

and At *= j ^ ^ (1-cosx) . 

^ atom 

The flux of electrons with velocities in the incremental 
range c^ , dc^ is 


c^n^f(c )dc 
e e ^ e’^ e 
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The rate at which these electrons pass through impact param- 
eter ring segments of area bdbd(j> centered about each heavy 
particle is given by Equation (C-1) : 


neHhf (Ce)CeO(x,Ce)sinCx)dxd<|)dc^ (C-1) 

The rate of energy loss due to collisions between electrons 

with velocities in c , dc and resulting in a deflection 

c c 

between x atid X'^’^ix is therefore, 

c^ 

J 3 — ^ Cl“COSx) 

enee e eoc <^atom 


Finally, the total energy loss rate due to collisions 
of all electrons deflected through all angles is 


2 e 

3 ^ 


atom 


( 

f2Tf r-TT 


I J o(x.Cg)sinx(l-cosx)dxd(j) 

j! 0 *^0 


dc 


The quantity in brackets is usually defined as the total 
elastic diffusion cross section CQ^j^) a^i^d is a function of 
Cg. Substituting the total cross section and expanding the 
differential vector velocity into its scalar spherical 
coordinate form yields the following expression for elec- 
tron atom Bremsstrahlung 
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f 2ir 
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Substituting the expression for the Maxwellian distribution 

unr 


, 7. 

and setting becomes: 


p — 2 e 




e^‘h ~TJY 
atom TT 


2kT^' 
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>"e 

V. / 





Since the integral is not significantly affected by 
varying the exponent of y by one (see Appendix D) , the 
exponent is reduced to 5 in order to satisfy Equation (D-5) 
for Jl = s = 1 , which is the gas kinetic collision integral 
for diffusion. The final form of the expression for elec- 
tron atom Bremss trahlung emission is then, 


128 


f ^2 'I 


fkT 1 


e h 

e 

1 c^d 

2m 

atom^ 

ej 


3/2 




where v 


1 

eh 




fkT^l 

® n 

e 

c*d . I 

2m 

atom-* 

e 

k. J 


V 


1 

eh 


is the collision frequency defined in Appendix D* 
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